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Growing and evolving: remarks for the 35" anniversary of the founding of 
Zoological Research 


The beginning of a new year is always a time for reflection, contemplation and optimism. For the past thirty-five years, 
since its founding in 1980, Zoological Research (ZR) has experienced many highs and lows, seeing great advances but 
also overcoming numerous hurdles. Zoological Research has made many valuable contributions to the progress of 
biological sciences in China and won prestige both at home and abroad. Currently, it serves as a key journal focusing on 
Genome Evolution and Genetic Diversity, Animal Ecology and Ethology, Primates and Animal Models of Human Diseases. 
The journal also publishes high quality papers related to taxonomy and DNA barcoding, developmental biology, physiology, 
biochemistry, immunology, and neuroscience. The growth achieved by ZH is not only due to our excellent readers and 
reviewers, but also to the evolution of our authors into renowned scholars. Accordingly, we would like to thank you all for 
your constant support and generous help. 

Zoological Research has always put the needs of its readers and authors first. To facilitate communication between our 
domestic colleagues and the international community, in 2013 to 2014 we successfully transitioned from a Chinese- 
language only journal, into a Chinese-English bilingual journal and finally into an internationally aimed English-language 
only journal. To help expand author view-points, we also included new Editorial and Letter to the Editor columns. 
Furthermore, the editorial board and editorial members continued to refine our publishing period, editing quality, typesetting, 
and journal marketing. In addition to ZR being selected as a high-quality publication for the dissemination of research 
findings, both the visibility and influence of ZR has continued to increase. 

As a respectable scientific publication in China, ZR had made impressive progress over the past few years, particularly 
in regards to its high academic value and leading role in advancing the quality of academic journals in China. In 2014, ZR 
was ranked among the top 300 “Outstanding S&T Journals of China", the second time ZR has been recognized by this 
prestigious award since 2008. In addition, four articles published in ZR between 2009 to 2013 were also included in the 
“Project of Frontrunner 5000” (F5000), which aims to promote scientific communication and the internationalization of 
Chinese S&T journals. 

The contributions from our colleagues have enabled our continued growth and advance. From a pool of more than 4600 
journals, ZR was among the top 10% of Chinese publications and was awarded the “The Highest / The International Impact 
Academic Journals of China” in Science, Technology and Engineering from 2012 to 2014 for its achievements in 
internationalizing life science research, as determined by its increasing international annual citation frequency and 
international impact factor. 

In September 2014, we were greatly honored to have a group of renowned academic experts from home and abroad join 
us as members of the editorial board. Thus, with dedicated assistance from our editorial team and first-rate reviewers, we 
will continue to expand and improve our publication, and in particular we welcome manuscripts on Genetics and Evolution, 
Primates and Animal Models of Human Diseases, and Biodiversity Conservation and Ecology. We also extend a special 
invitation to our colleagues to serve as single-issue Guest Editors, and if any author would like to advance a topic that is 


prospective, controversial and worth exploring, please feel free to contact us. 
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Once again, we would like to thank you all for your enduring support and faith. Our continued growth and evolution 


together will ensure that ZR remains a respected publication platform that provides ever greater possibilities to you. 


Sincerely yours, 


~ 


Yong-Gang YAO, Editor-in-Chief 


Kunming Institute of Zoology, Chinese Academy of Sciences, Kunming 650223, China 
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Yun ZHANG, Executive Editor-in-Chief 


Kunming Institute of Zoology, Chinese Academy of Sciences, Kunming 650223, China 
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ABSTRACT 


Ruminant stomach lysozyme is a long established 
model of adaptive gene evolution. Evolution of 
stomach lysozyme function required changes in the 
site of expression of the lysozyme c gene and 
changes in the enzymatic properties of the enzyme. 
In ruminant mammals, these changes were 
associated with a change in the size of the 
lysozyme c gene family. The recent release of near 
complete genome sequences from several 
ruminant species allows a more complete 
examination of the evolution and diversification of 
the lysozyme c gene family. Here we characterize 
the size of the lysozyme c gene family in extant 
ruminants and demonstrate that their pecoran 
ruminant ancestor had a family of at least 10 
lysozyme c genes, which included at least two 
pseudogenes. Evolutionary analysis of the ruminant 
lysozyme c gene sequences demonstrate that each 
of the four exons of the lysozyme c gene has a 
unique evolutionary history, indicating that they 
participated independently in concerted evolution. 
These analyses also show that episodic changes in 
the evolutionary constraints on the protein 
sequences occurred, with lysozyme c genes 
expressed in the abomasum of the stomach of 
extant ruminant species showing the greatest levels 
of selective constraints. 


Keywords: Lysozyme c; Ruminants; Gene family; 
Gene duplication; Concerted evolution; Mosaic 
evolution 


INTRODUCTION 


Ruminant mammals such as cow, sheep, and deer, rely on 
foregut fermentation to extract nutrients from their diet of 
plant material (Clauss et al, 2010; Janis, 1976; Mackie, 2002; 
Stevens & Hume, 1998). Foregut fermentation, by bacteria 
and microbes, produces short chain fatty acids that are 
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absorbed through the stomach wall and provide energy for 
the ruminant animals; however, the microbial population 
responsible for this fermentation incorporates many of the 
other nutrients, such as nitrogen based compounds, into their 
own growing populations (Mackie, 2002; Stevens & Hume, 
1998). To extract these essential nutrients from the microbial 
population, ruminant animals must break open these bacterial 
and microbial cells, to release their contents, to allow the 
stomach digestive enzymes in the abomasum to extract 
nutrients from their contents (Stevens & Hume, 1998). Since 
bacterial cells are typically resistant to mammalian digestive 
enzymes, ruminant species have recruited the anti-bacterial 
enzyme, lysozyme c, to break open these cells (Callewaert & 
Michiels, 2010; Dobson et al, 1984; Irwin et al, 1992; Mackie, 
2002; Prager & Jollés, 1996). Recruitment of lysozyme c as a 
digestive enzyme has occurred at least twice within mammals, 
on the lineages leading to the ruminant artiodactyls and the 
leaf-eating monkeys (Dobson et al, 1984; Stewart et al, 1987; 
Stewart & Wilson, 1987), with a similar recruitment of a calcium- 
binding lysozyme occurring in the hoatzin, a leaf-eating bird 
(Kornegay et al, 1994; Kornegay, 1996). 

Recruitment of lysozyme c to become a digestive enzyme 
required changes both in the site of expression of the gene 
encoding this enzyme and in the amino acid sequence of the 
enzyme to allow function in the acidic stomach (Dobson et al, 
1984; Irwin et al, 1992; Irwin, 1996; Prager, 1996). The 
major site of expression of lysozyme in mammals is 
macrophages, but it also secreted into some body fluids 
(such as tears), where it participates in host defense against 
bacterial infection (Callewaert & Michiels, 2010; Prager & 
Jollés, 1996; Short et al, 1996). The molecular basis for the 
recruitment of expression, at high levels, of lysozyme c in 
stomach cells is unknown. Typical mammalian lysozyme c 
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enzymes function in an environment at a neutral pH, and 
one that is free of digestive enzymes (Callewaert & Michiels, 
2010; Prager & Jollés, 1996; Prager, 1996). Lysozyme c 
function in the abomasum of the stomach of ruminant 
animals, to digest bacterial cell walls, required adapting the 
lysozyme c protein sequence to function at an acidic pH and 
becoming resistant to the actions of stomach digestive 
enzymes and acids found in the abomasum (Dobson et al, 
1984; Jollés et al, 1989; Prager, 1996). A number of 
convergent amino acid changes were seen between the 
lysozyme c sequences that have adapted for function in the 
stomachs of the langur, a leaf-eating monkey, and ruminants, 
have been identified and presumed to account for much of 
the functional adaptation (Stewart & Wilson, 1987; Stewart 
et al, 1987; Swanson et al, 1991; Prager, 1996). Some of 
these adaptive changes include replacement of lysine 
residues with arginine, which removes potential cleavage 
sites for digestive enzymes found in the stomach, and the 
loss of an aspartate-proline dipeptide, which is an acid-labile 
peptide bond (Jollés et al, 1989; Prager, 1996; Stewart & 
Wilson, 1987; Stewart et al, 1987; Swanson et al, 1991). 
These putative adaptive amino acid replacements are 
inferred to occur early in ruminant evolution, and thus may 
parallel the origin and evolution of the ruminant lifestyle 
(Irwin et al, 1992; Irwin, 1996). 

Recruitment of lysozyme c to a digestive role in ruminants 
is associated with an expansion of the size of the lysozyme 
c gene family (Jiang et al, 2014; Irwin & Wilson, 1989; Irwin 
et al, 1989, 1992). Most mammals have only one or a few 
lysozyme c genes, while ruminant species have 10 or more 
(Callewaert & Michiels, 2010; Irwin & Wilson, 1989; Irwin et 
al, 1989, 1996; Prager & Jollés, 1996; Irwin et al, 2011; 
Jiang et al, 2014). The lysozyme c gene family of the cow 
has been better characterized than those of most other 
ruminants, where it was found that only some of the genes 
are expressed in the abomasum, while others retain more 
ancestral type of roles (Irwin & Wilson, 1989; Irwin et al, 
1993; Irwin, 2004). Similar observations have been made for 
the lysozyme c genes of sheep (Jiang et al, 2014). Several 
lysozyme c proteins, and their cDNAs, have been 
characterized from the abomasums of the cow, sheep and 
deer (Dobson et al, 1984; Jollés et al, 1989; Irwin & Wilson, 
1989, 1990). Intriguingly, phylogenetic analysis of the coding 
and 3' untranslated portions of the lysozyme c cDNA 
sequences yielded different trees, with the coding 
sequences implying duplications of the genes on each 
species lineage and the 3 untranslated region indicating 
more ancient duplications before the divergence of these 
species (Irwin & Wilson, 1990). Selection at the protein level 
(e.g., lineage-specific adaptation of the protein sequences) 
does not explain the differences in the phylogenies for the 
two regions, as synonymous difference (those that do not 
change the coding potential) also yield the same 
conclusions. It was concluded that the differences in the 
phylogenies was due to concerted evolution, mediated by 
gene conversion, acting on the coding sequences, while the 
3' untranslated regions only experienced divergent evolution 
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(Irwin & Wilson, 1990; Irwin et al, 1992; Irwin, 1996; Wen & 
Irwin, 1999; Yu & Irwin, 1996). Characterization of the 
genomic sequences of lysozyme genes expressed in the 
abomasum from the cow and sheep suggested that the 
concerted evolution was limited to only the coding exons, 
and did not involve the intronic sequences separating these 
exons (Irwin et al, 1993; Wen & Irwin, 1999). An analysis of 
larger number of lysozyme mRNA sequences, including 
genes that are not expressed in the stomach, suggested that 
some of the genes expressed in non-stomach tissues might 
have also experienced concerted evolution (Irwin, 1995, 
2004; Takeuchi et al, 1993). 

The previous analyses were largely limited to lysozyme c 
genes expressed in ruminant species. With the recent 
completion of draft genomic sequences from several 
ruminant species (including, cow, yak, zebu, goat, Tibetan 
antelope, and sheep: Canavez et al, 2012; Dong et al, 2013; 
Ge et al, 2013; Jiang et al, 2014; Qiu et al, 2012; Zimin et al, 
2009) it now possible to more completely characterize the 
complete complement of lysozyme c genes (including 
genenes that are not expressed) in the genomes of these 
species and examine the molecular evolution of these genes. 
Here we describe the lysozyme c gene complements of the 
cow and several other ruminant species. The lysozyme c 
gene cluster has largely been maintained within true 
ruminant (Infraorder Pecora) species. Analysis of these 
sequences shows that the ancestor of cow, sheep, and 
goats had 10 lysozyme c genes, several of which were 
pseudogenes that were retained by diverse species. The 
exons of the lysozyme c genes have differing evolutionary 
histories, suggesting that concerted evolution acted 
independently on each exon. 


MATERIALS AND METHODS 


Database searches 

Previous searches of mammalian genomes indicated that 
the cow genome had about 12 lysozyme c genes located in 
a cluster on cow chromosome 5, many of which were 
incompletely annotated in the Ensembl assembly (Irwin et al, 
2011). To better characterize the lysozyme c gene cluster in 
the cow (Bos taurus) genome we used the Blast algorithm 
(Altschul et al, 1990) to search the UMD 3.1 cow genome 
assembly (from Ensembl release 75 in June 2014; 
http:/Awww.ensembl. org/index.html) with known and predicted 
cow lysozyme c cDNA and protein sequences. Lysozyme c 
genes from the sheep (Ovis aries; Oar_v3.1), pig (Sus 
scrofa; Sscrofa10.2), bottlenose dolphin (Tursiops truncates; 
Turtru1), dog (Canis lupus familiaris; CanFam3.1), panda 
(Ailuropoda melanoleuca; AilMel1), horse (Equus caballus; 
EquCab2), and rhinoceros (Ceratotherium simum simum; 
CerSimSim1 - preEnsembl) genomes from the Ensembl 
database were characterized by the approaches described 
above. A similar search strategy was used to identify 
lysozyme c genes in the yak (Bos grunniens), zebu (Bos 
indicus), water buffalo (Babalus babalis), Tibetan antelope 
(chiru; Pantholopus hodgsonii), goat (Capra hircus), aplaca 


(Vicugna pacos), minke whale (Balaenoptera acutorostrara 
scammoni), killer whale (Orcinus orca), Yangtze River 
dolphin (Lipotes vexillifer), and sperm whale (Physeter 
catodon) genomes from the NCBI Genomes (chromosome), 
Whole-genome shotgun contigs (wgs), and Nucleotide 
collection (nr/nt) databases (http://blast.ncbi.nlm.nih.gov/Blast. 


cgi). 


Genomic alignments and assignment of orthology 
Genomic sequences encompassing lysozyme c genes were 
downloaded from the Ensembl and NCBI databases. Intron- 
exon boundaries of the 4 exons and the 5' and 3' flanking 
sequences of the new lysozyme c genes were annotated 
based on genomic alignments of genes using MultiPipMaker 
(Schwartz et al, 2000, 2003), using previously characterized 
artiodactyl lysozyme c genes (Irwin et al, 1993, 1996; Irwin, 
1995; Yu & Irwin, 1996; Wen & Irwin, 1999) as guides. Gene 
neighborhood organization was assessed as previously 
described for lysozyme c genes (Irwin et al, 2011) with the 
flanking Yeats4 and Cpsf6 genes identified using Blast. 
Ruminant lysozyme c genes were named based on 
orthology (based on phylogeny, see below) and genomic 
location. Genes present in the common ancestor of sheep 
and cow were numbered (Lyz1-Lyz10), while lineage- 
specific duplicates have a letter (a-c) that follows the gene 
number. The alpaca lysozyme c genes were numbered 
arbitrarily, thus their numbers do not indicate orthology with 
the ruminant genes. All other species examined here have a 
single copy lysozyme c gene. 


Phylogenetic analysis 

Predicted protein coding sequences for lysozyme c cDNA 
sequences, extracted from the genomic alignments, were 
aligned with Muscle (Edgar, 2004) as implemented in 
Mega6.06 (Tamura et al, 2013). Alignments were edited 
manually to insert gaps to maintain open reading frames 
(due to the presence of frame shifting insertions in some 
pseudogenes). Phylogenetic trees were constructed by the 
maximum likelihood, neighbor-joining and parsimony 
methods using Mega6.06 (Tamara et al, 2013). Alternative 
phylogenetic hypothesis, derived from the phylogenies of 
the different exons, were tested using Tree-puzzle (Strimmer 
& von Haeseler, 1996) as implemented on the Mobyle 
(QPasteur web site (http://mobyle. pasteur.fr/cgi-bin/portal.py? 
welcome; Néron et al, 2009). 


RESULTS AND DISCUSSION 


Number and organization of lysozyme c genes in the cow 
genome 

Analyses of genomic Southern blots had concluded that 
there were about 10 lysozyme c genes in the cow genome 
(Irwin & Wilson, 1989; Irwin et al, 1989), with many of these 
genes clustered on chromosome 5 (Gallagher et al, 1993). A 
recent (2011) search of the Btau4.0 (2™ release, assembled 
2007) of the cow genome sequence assembly identified 12 
lysozyme c genes on chromosome 5 of the cow genome 


(Irwin et al, 2011). The genes identified in this search 
account for all of the previously characterized cow lysozyme 
c cDNA and protein sequences (Irwin et al, 2011). To better 
characterize the cow lysozyme c gene cluster, we searched 
the most current version (UMD3.1 — 3" release) of the cow 
genome assembly (assembled 2009; Zimin et al, 2009) with 
Blast using the previously characterized lysozyme c cDNA 
and protein sequences. Our new searches identified a total 
of 14 lysosome c genes, 11 of which were annotated by 
Ensembl as genes (Table 1 and Figures 1, 2, and S1, 
supporting information at http://www.zoores.ac.cn/). The 
difference in the number of intact genes identified by the 
searches of the two different genome assemblies (12 in 
Btau4.0 and 11 in UMD3.1) is due to the earlier Btau4.0 
assembly containing two copies of the tracheal lysozyme c 
gene (CowC and CowD in Irwin et al, 2011) while the most 
current assembly UMD 3.1 contains only single copy of this 
gene (here named Lyz2b). 

In addition to the 11 annotated genes, each of which is 
composed of 4 exons consistent with the structure of a 
typical mammalian lysozyme c gene (Irwin et al, 1996; 
Callewaert & Michiels, 2010), Blast hits were found to map 
to additional locations that were distant from the annotated 
genes. Examination of these Blast hits suggested that they 
belong to three partial genes, which had not previously been 
annotated, with each being composed of only two, not four, 
exons (Table 1 and Figures 1, 2, and S1). The newly 
identified partial Lyz3a and Lyz3c genes contain exons 1 
and 2 and exons 3 and 4, respectively, but are separated 
from each other by the Lyz3b gene that contains all 4 coding 
exons (Table 1 and Figures 1 and 2). No sequences similar 
to the missing exons were found near the Lyz3a and Lyz3c 
genes. The third partial gene, Lyz9, contains only exons 1 
and 4, with the sequences between these exons showing no 
similarity to exons 2 or 3 of other lysozyme c genes. Most of 
the lysozyme c genes have the same orientation (annotated 
as the minus strand), but 4 of the 14 are on the opposite 
strand, indicating that the origin of this gene family is not just 
a simple series of tandem gene duplications. The cow 
genome, therefore, was found to contain 14 identifiable 
lysozyme c genes (Table 1 and Figures 1,2, and S1). 

Pairwise sequence comparisons revealed that the DNA 
sequence identities of the coding sequences among the 14 
genes ranged from 74.896 to 97.596, with most pairs 
showing 8096-9096 identity (Table 2). The similarity between 
Lyz3a and Lyz3c could not be measured, as there is no 
overlap between these two genes (Lyz3a has exons 3 and 4, 
while Lyz3c has exons 1 and 2, see Table 1). These two 
partial genes were most similar to the Lyz3b gene, showing 
greater than 96% identify in the coding sequence (Table 2), 
raising the possibility that they are recent gene duplicates. 
The Lyz3a, Lyz3b, and Lyz3c genes are also adjacent to 
each other in the genome, suggesting that the Lyz3a and 
Lyz3b gene were generated by partial tandem duplications 
of different parts of the Lyz3b gene (Figure 1). The partial 
gene Lyz9 did not show particularly strong similarity to any 
other specific cow lysozyme c gene (Table 2), suggesting 
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that it is not a product of a very recent segmental duplication 
event (Liu et al, 2009; Seo et al, 2013). Among the intact 
lysozyme genes, the coding sequence of the genes 
encoding the lysozymes expressed in the abomasum (Irwin 


& Wilson, 1989; Irwin et al, 1993) Lyz5/Lyz6/Lyz7 share 
about 97% identity and the Lyz2a/Lyz2b/Lyz2c genes, which 
includes the tracheal lysozyme gene (Takeuchi et al, 1993), 
share about 96% identity (Table 2). The high level of identity 








Table 1 Locations of cow lysozyme c genes 
Oth Sites of 
Gene Chromosome Strand Bases Ensembl gene ID e Unigene bae ! Functional 
name expression 
; 44 421 190- 7 
Lyz1 5 Minus 44 426 117 ENSBTAG00000011941 Milk Bt.67194: 448 ESTs Rumen, omasum Intact 
Lyz2a 5 Plus ppc uA ENSBTAGO00000022971 None? Intact 
y 44 426 117 
44 443 587- 
Lyz2b 5 Plus ENSBTAGO00000000198  Trachael X Bt.64327: 121 ESTs Rumen Intact 
44 448 198 
44 489 738- 
Lyz3a 5 Minus NA (exons 3 and 4)? None? Pseudogene 
44 491 389 
44 502 011- 
Lyz3b 5 Minus ENSBTAG00000039170 wNS4 None? Pseudogene 
44 507 108 
44 521 088- 
Lyz3c 5 Minus NA (exons 1 and 2)? None? Pseudogene 
44 522 921 
44 533 912- : . 
Lyz2c 5 Plus ENSBTAGO00000020564 Bt.105675: 31 ESTs Rumen, intestine Intact 
44 538 847 
44 554 351- 
Lyz4 5 Plus 44 559 113 ENSBTAG00000026323 Intestinal — Bt.49176: 162 ESTs Intestine Intact 
: 44 573 344- 
Lyz5 5 Minus ENSBTAG00000026088 Stomach2 Bt.29367: 363 ESTs Abomasum Intact 
44 578 495 
: 44 599 815- 
Lyzé6 5 Minus 44 607 109 ENSBTAGO00000046511 Stomach 1 Bt.89770: 102 ESTs Abomasum Intact 
: 44 631 778- 
Lyz7 5 Minus 44 638 267 ENSBTAG00000046628 Stomach3 Bt.80498:74 ESTs Abomasum Intact 
] 44 652 535- b 
Lyz8 5 Minus 44 656 613 ENSBTAG00000026322 None Pseudogene 
; 44 673 817- : b 
Lyz9 5 Minus 44 676 391 NA (exons 1 and 4) None Pseudogene 
Lyz10 5 Minus 74713299- ^ ENSBTAG00000026779 Kidney ^ BL64645:61ESTs MPhoreticular rect 
14 44 720 946 d binds blood 
? — Not annotated as a gene by Ensembl. 
? — No ESTs with greater than 9596 sequence identity identified in the NCBI EST database. 
Lyz2a Lyz2b 73220 Lyz4 
m m m m 
NMM W VW Wo M WM Ww Wy IMI MMMA 
Yeats4 Lyz1 Lyz3b Lyz5 27026 2027 Lyz9 Lyz10 Cpsf6 
Lyz3a 70230 7028 
44.3 Mb 44.4 Mb 44.5 Mb 44.6 Mb 44.7 Mb 44.8 Mb 
l ㅣ l ㅣ | ㅣ ㅣ ㅣ | 
500 kb > 








Chromosome 5 


Figure 1 Organization of lysozyme c genes in the cow genome 

Schematic of the arrangement of lysozyme c genes. And their neighbors, in the cow genome. Vertical lines represent exons, with splicing indicated by the 
lines joining the exons. Gene names are indicated above (plus strand) or below (minus strand) indicating strand with coding potential. Sizes of genes and 
distances are proportional. The genes are located between 44.3 and 44.8 Mb on chromosome 5. 
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Lyzl MKALLIVGLL LLSVAVQG KKFQRCELAR TLKKLGLDGY RGVS-LANWV CLARWESNYN TRATNYNRGD KSTDYGIFQI 
Lyz2a 365 3 K...-....M ...KG..G.. 
Lyz2b Koominen M cos EG. Ss. 
Lyz2c K = DM .«.KG..8.. 
Lyz3a | ---------- -------- ---------- ---------- ---------- ---------- ------- 
Lyz3b 00020 Ras: aS BRE 
Lyz3c K.I.-..K.M ...8..R... 
Lyz4 K..2292. aM aa TYG Ray 
Lyzs Cea area) mers 248 = 
Lyz6 Ecosse ba TR o Gaa 
Lyz7 E E es ees 
Lyz8 Koo ROL LITERIS 080. 
Lyz9 -...?- ---------- ---------- ---------- 
Lyz10 I.- Ms ODE Eo A [ope 
$ # # 
100 129 
Lyzl NSRWWCNDGK TPKAVNACRI PCSALLKDDI TOAVACAKRV VRDPOGIKAW VAWRNKCONR DLRSYVOGCRV 
Lyz2a uu "CEP P KI .SQ-L.BET.. ace Kis eae es. s EEEE A 
18020: ^ 30698 xe x GV Bue ieedee eu eode KI i809... oxexcBRer oc. Tes Kasas 
Lyz2c Bes vel ete as GM. Sia ra been ees Cav TR D KI 490-2 D. ara Ke Ni Beas: se 인게 we G. 
Lyz3a ---------- ---------- ----- Quot eS B IK.DIdES..2 R2 SER Ade edt Or 2 eRe ons. G. 
Lyz3b 고도. aa 시가: AONNE + dyes) aie B.oK. S.r VR 원. [^ ~ D.G 
Lyz3c PE eres «BB nea 시 "or 
Lyz4 SUR eset a, Maud gus GGV SIL I 내에 새우 e ag Sue be TI .SR-...T K...R... VS..IR..KL 
Lyz5 LU c ..N..DG.HV S..E.MEN.. AK...... HI .SE- "ID. sc eo KSHARDE: 2MS..zE.TD 
Lyz6 PEE ata Vs ..N..DG.HV S..E.MEN.. AK...... QI .SE-...T.o. aea KSH- RDH .VS-.lE..TD 
Lyz7 EI SEP ..N..DG.HV S..E.MEN.. AK...... HI .SE- T das RSILGRIEH- Veces TL 
Lyz8 vel ade um der DG.PV SH.K.MGN.. AK...... 1 SE- T.. ...KSH.RDH .VS...E..TL 
Lyz9 ---------- ---------- c--------- ---------- --------- P^ o sRNOMSIQ. acosan Diss 
Lyz10 oe Hed ee tans ves Hb 2dRborsQkee maestru iti Serisi Ride ea SH as0 ai Des Lewe Ss 
Figure 2 Amino acid sequences predicted by cow lysozyme c genes 
Sequences of predicted lysozyme c proteins from the cow genome are shown in single letter code, with differences from the Cow Lyz1 sequence shown 





and identities indicated by dots (.). Sequence is numbered avove the sequences from the N-terminus of the Lyz1 sequence, with the signal peptide 
numbered backwards and in italics. Dashes (-) indicate gaps introduced to maximize alignment. Question marks (?) indicate incomplete codons due to 
missing sequence. Residues involved in disulfide bridging ($) and active site residues (#) are marked below the sequences. Residues marked in red are 
likely damaging pseudogenes and disrupt initiation, disrupt disulfide bridging or introduce stop codons. Asterisks (*) indicate inframe stop codons. Xs 
refer to codons that have less than 3 bases and thus cause frame shifts. The initiation codon of Lyz8 is not methionine (M). 


Table2 Pairwise percent DNA sequence identity between cow lysozyme c coding sequences 
Lyz2a Lyz2b Lyz3a Lyz3b Lyz3c Lyz2c Lyz4 Lyz5 Lyz6 Lyz7 Lyz8 Lyz9 Lyz10 
Lyz1 89.4 88.7 91.3 91.3 90.7 89.0 86.5 82.2 83.1 82.9 80.9 85.6 89.0 





Lyz2a 96.2 87.4 85.3 84.4 96.6 90.1 84.7 85.6 85.4 82.8 84.7 84.5 
Lyz2b 88.1 85.3 84.1 96.8 90.5 84.9 85.8 85.4 82.8 83.2 84.5 
Lyz3a 96.4 NA? 88.1 83.0 74.8 75.6 76.3 76.3 77.3 89.9 
Lyz3b 96.3 85.3 85.3 79.0 79.5 79.7 77.6 84.2 86.5 
Lyz3c 84.1 86.0 82.1 82.4 82.4 79.5 86.0 86.4 
Lyz2c 89.6 84.9 85.8 85.6 82.8 83.7 84.7 
Lyz4 84.7 84.7 85.4 83.0 83.7 83.3 
Lyz5 96.4 97.5 92.2 83.2 82.9 
Lyz6 96.8 92.2 84.2 82.0 
Lyz7 92.4 84.7 82.9 
Lyz8 82.7 80.0 
Lyz9 85.6 


? — These two genes do not overlap. 


shared by these sets of genes suggests that these triplets Lyz5/Lyz6/Lyz7 are adjacent to each other and are in the 
are products of recent segmental duplication / gene same orientation (Table 1; Figure 1), thus could be 
duplication events (Liu et al, 2009; Seo et al, 2013). generated by a simple series of tandem gene duplication 
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events. A more complicated duplication history is needed to 
explain the diversification of the Lyz2a/Lyz2b/Lyz2c genes. 
While Lyz2a and Lyz2b are in tandem, several other 
lysozyme c genes (Lyz3a, Lyz3b, and Lyz3c) are located 
between the Lyz2a/2b gene pair and the Lyz2c gene (Table 
1 and Figure 1). 

Since the three partial genes Lyz3a, Lyz3b, and Lyz9 do 
not contain all four coding exons; they cannot predict intact 
open reading frames. In addition to the missing exon 
sequence, all three of these genes also contain in frame 
stop codons or frameshifts that also would prevent 
translation (Figure 2). Among the lysozyme c genes 
possessing all four exons, two, Lyz3b and Lyz8, fail to 
predict intact open reading frames (Figure 2). Lyz3b, which 
was previously called lysozyme wNS4 (Irwin, 1995) contains 
a frameshift, which is shared with Lyz3a, which prevents 
translation of the reading frame, while Lyz8 contains both in 
frame stop codons and a replacement at the initiating codon 
(Figure 2). Thus of the 14 lysozyme c genes found in the 
cow genome, only 9 potentially encode functional lysozyme 
c proteins. To further investigate the functional potential of 
these lysozyme c genes we searched for evidence of expre- 
ssion for all 14 cow lysozyme c genes in the NCBI expre- 
ssed sequence tag (EST) database. ESTs were found for 
only 8 of the 9 intact genes, and for none of the 5 
pseudogenes (Table 1). While Lyz2a has an intact open 
reading frame (Figure 2), no ESTs highly similar (>98% 
identity) to it were found in the NCBI database (Table 1), 
raising the possibility that this gene is not expressed. 

Many of the cow lysozyme c gene annotations in the 
Ensembl database do not include 5' and/or 3' untranslated 
sequences. Since previous work had shown that the 5' and 
3' untranslated sequences of known lysozyme c genes from 
diverse mammalian species have considerable sequence 
similarity (Irwin & Wilson, 1989; Irwin, 1995, 2004), we used 
this similarity to predict the extent of these regions for each 
gene (see Figure S1) from alignments generated by 
MultiPipMaker (Schwartz et al, 2000, 2003). Complete 5' 
untranslated regions could be predicted for all of the 
lysozyme c genes that had exon 1, however the full 3' 
untranslated regions could not be predicted for all exon 4 
sequences, as the 3' end of the 3' untranslated region could 
not be found for the cow Lyz3a and Lyz3b genes (see 
Figure S1). This observation is consistent with an earlier 
failure to identify homologous sequences for the entire 3' 
untranslated region of the cow wNS4 (Lyz3b) gene (Irwin, 
1995, 2004). 


Lysozyme c genes in other ruminant genomes 

To better characterize the evolutionary history of the 
ruminant lysozyme c genes, we identified lysozyme c genes 
in the genomes of other ruminant species and their close 
relatives (Table 3 and Figures 3 and S1). As expected, from 
previous work (Callewaert & Michiels, 2010; Irwin et al, 1989, 
1996, 2011; Prager & Jollés, 1996), only a single lysozyme c 
gene was found in the genomes of carnivores (dog, Canis 
lupus familiaris; and panda, Ailuropoda melanoleuca) and 
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perrisodactyls (horse, Equus caballus, and rhinoceros, 
Ceratotherium simum simum) (Table 3 and Figures 3 and 
S1). The single lysozyme c gene in the outgroup species is 
located between the Yeats4 and Cpsf6 genes (Figure 3), as 
found in most other mammalian species (Irwin et al, 2011). 
This ancestral mammalian genomic arrangement has been 
retained in the cow, with the amplification of the lysozyme c 
genes occurring between the Yeats4 and Cpsf6 genes (Irwin 
et al, 2011) (Figures 1 and 3). The tylopod lineage (e.g., 
camels and alpacas) represents one branch of the earliest 
divergence within artiodactyls (Morgan et al, 2013; 
Romiguier et al, 2013), with these species being 
pseudoruminants with a simpler multi-chambered stomach 
than the true ruminants (Clauss et al, 2010; Janis, 1976; 
Mackie, 2002). Searches of the alpaca (Vicugna pacos) 
genome in the Ensembl database identified three genomic 
sequences encoding partial lysozyme c gene sequences, 
indicating that multiple lysozyme c genes exist in this 
genome (results not shown). Searches of the NCBI 
Genomes (chromosomes) database identified an updated 
larger genomic contig that predicted 4 complete lysozyme c 
genes (and included all of the gene sequences found in the 
Ensembl alpaca genome assembly) at one end of a contig 
sequence (Table 3 and Figure 3). The Yeats4 gene was 
found to be adjacent to one side of the lysozyme c gene 
cluster, however no genes were found on the other side of 
the lysozyme c gene cluster in this genomic contig (Figure 
3). The presence of the Yeats4 gene adjacent to the alpaca 
lysozyme c genes suggests that a similar genomic 
neighborhood exists in alpaca, but since the lysozyme c 
genes were at one end of the genomic contig it is possible 
that additional unsequenced lysozyme c genes may exist in 
the alpaca genome. The pig (Sus scrofa) is a representative 
of the family Suidea, which is the next diverging lineage 
within artiodactyls (Morgan et al, 2013; Romiguier et al, 
2013). As expected, and previously reported (Irwin et al, 
1989; Yu & Irwin, 1996), only a single lysozyme gene is 
found in this species (Table 3, Figures 3 and S1). As 
previously reported (Irwin et al, 2011), the genomic 
neighborhood surrounding the pig lysozyme gene differs 
from that of other mammals, raising the possibility that this 
genomic area has experienced recombination (Figure 3). 
Cetaceans (e.g., whales and dolphins) fall within artiodactyls, 
thus yielding cetartiodactyla (Morgan et al, 2013; Romiguier 
et al, 2013). A single lysozyme c gene was identified in all 
five cetacean (bottlenose dolphin, Tursiops truncates; minke 
whale, Balenoptera acutorostrata scammoni; killer whale, 
Orinus orca; Yangtze river dolphin, Lipotes vexilllifer, and 
sperm whale, Physeter catodon) genomes (Table 3 and 
Figure S1), which is found in genomic location consistent 
with the ancestral genomic organization (Figure 3). 

Pecoran artiodactyls (cow, sheep, deer, and relatives) are 
true ruminants with a stomach composed of four chambers 
(Clauss et al, 2010; Janis, 1976; Mackie, 2002). In addition 
to the sheep (Ovis aires) genome (Jiang et al, 2014), which 
is available from Ensembl, genome sequences of 5 other 
pecoran ruminant species (yak, Bos grunniens (Qiu et al, 


Table 3 Locations of lysozyme c genes in diverse artiodactyls and relatives 


Chromosome / 


Ensembl gene ID / 





Gene Strand Bases d Missing exons 
scaffold NCBI accession 
Yak (Bos grunniens) 
Lyz1 NW. 005394307 Minus 5 211-9 033 XM 005901148 
Lyz2 NW. 005394307 Plus 75 876-81 020 NA? 1 
Lyz3 NW. 005394198 Plus 7 980-12 653 XM 005900299 
Lyz4 NW. 005394198 Minus 27 041-32 422 XM. 005900300 
Lyz5 NW. 005394198 Minus 53 712-61 005 XM. 005900301 
Lyz6 NW. 005394198 Minus 83 939-90 419 XM. 005900302 
Lyz7 NW_005394198 Minus 104 585-108 706 NA? 
Lyz8 NW. 005392857 Minus 16 134-18 709 NA? 2,3 
Lyz9 NW_005392857 Minus 61 437-66 594 XM_005886999 
Lyz10 NW_005394307 Minus 5 211-9 033 XM 005901148 
Zebu (Bos indicus) 
Lyz1 AGFL01046860 Minus 920-2 950 NA? 3, 4° 
Lyz2a AGFL01046876 Plus 10 835-14 065 NA? 
Lyz2b AGFL01046877 Plus 11 353-15 933 NA? 
Lyz3 AGFL01046880 Minus 12 055-17 431 NA? 
Lyz2c AGFL01046883 Plus 3 743-8 620 NA? 
Lyz4 AGFL01046890 Plus 880-5 642 NA? 
Lyz5 AGFL01046890 Minus 19 872-25 024 NA? 
Lyz6 AGFL01046892 Minus 13 424-20 718 NA? 
Lyz7 AGFL01046895 Minus 668-7 157 NA? 
Lyz8 AGFL01046896 Minus 13 855-17 933 NA? 
Lyz9 AGFL01046897 Minus 889-1 024 NA? 2,3,4° 
Lyz10 AGFL01046900 Minus 2 248-9 934 NA? 
AGFL01046860 Minus 920-2 950 NA? 3, 4° 
Water buffalo (Bubalus bubalis) 
Lyz1 NW_005784949 Minus 16 204-27 136 XM_006058377 
Lyz2 NW_005785126 Plus 13 244-16 121 NA? 4° 
Lyz3 NW_005785126 Minus 148 364-153 473 NA? 
Lyz4 NW_005785126 Plus 173 780-178 474 XM_006064264 
Lyz5 NW. 005785126 Minus 194 055-199 224 XM 006064265 
Lyz6 NW. 005785126 Minus 220 532-227 898 XM 006064266 
Lyz7 NW_005785126 Minus 245 385-151 911 XM_006064267 
Lyz8 NW. 005785126 Minus 269 603-271 738 NA? 
Lyz9 NW. 005785126 Minus 283 880-286 454 NA? 2,3 
Lyz10 NW. 005785126 Minus 323 066-328 240 XM 006064268 
NW 005784949 Minus 16 204-27 136 XM 006058377 
Tibetan antelope / Chiru 
(Pantholops hodgsonii) 
Lyz1 NW. 005806187 Minus 948 632-954 099 XM 005957102 
Lyz2 NW. 005806187 Plus 1 005 582-1 010 034 XM 005957103 
Lyz3 NW. 005806187 Minus 1 045 148-1 051 559 XM 005957104 
Lyz4 NW. 005806187 Plus 1 074 291-1 079 230 NA? 2, 3? 
Lyz5 NW_005811703 Minus 9 411-10 957 XM_005968144 3, 4° 
Lyz6 NW. 005811703 Minus 31 570-36 493 XM 005968125 
Lyz7 NW_005811703 Minus 64 104-69 686 XR_318952 
Lyz8 NW. 005811703 Minus 83 106-87 330 NA? 
Lyz9 NW_005811703 Minus 102 887-105 253 NA? 2,3 
Lyz10 NW_005811703 Minus 168 707-174 027 XM_005968126 
Goat (Capra hircus) 
Lyz1 NW_005100667 Minus 6 548 409-6 554 516 XM_005680191 
Lyz2 NW_005100667 Plus 6 607 855-6 612 306 XM_005680189 
Lyz3 NW_005100667 Minus 6 638 119-6 644 658 NA? 
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Continued 





Guns Chromosome / Strandi - ‘Bases Ensembl gene ID/ Missing 
scaffold NCBI accession exons 

Lyz4 NW_005100667 Plus 6 670 972-6 67 2304 NA? 1,2° 

Lyz5 NW. 005100667 Minus 6 689 705-6 694 702 XM 005680192 

Lyz6 NW. 005100667 Minus 6 715 083-6 719 993 NM, 001287566 

Lyz7 NW_005100667 Minus 6 741 973-6 747 565 NA? 

Lyz8 NW. 005100667 Minus 6 761 023-6 765 243 XM 005680235 

Lyz9 NW. 005100667 Minus 6 783 196-6 785 752 NA? 2,3 

Lyz10 NW. 005100667 Minus 6 834 500-6 839 640 NM 001285711 

Sheep (Ovis aries) 

Lyz1 3 Minus 150 165 176-150 170 352 ENSOARGO00000020393 

Lyz2 3 Plus 150 225 205-150 229 630 ENSOARGO00000020417 

Lyz3 JH921983.1 Minus 4 380-1 032 ENSOARGO00000000543 

Lyz4 3 Plus 150 266 228-150 270 937 ENSOARGO00000020429 

Lyz5 3 Minus 150 288 480-150 293 529 ENSOARGO00000020393 

Lyz6 3 Minus 150 313 875-150 318 810 ENSOARGO00000020439 

Lyz7 3 Minus 150 342 914-150 348 498 NA? 

Lyz8 3 Minus 150 362 122-150 366 351 ENSOARGO00000020476 

Lyz9 3 Minus 150 385 062-150 387 578 NA? 2,3 

Lyz10 3 Minus 150 434 372-150 439 510 ENSOARGO00000020515 

Pig (Sus scrofa) 

Lyz 5 Minus 36 179 949-36 185 488 ENSSSCG00000000492 5 

Alpaca (Vicugna pacos) 

Lyz1 NT_167289.2 Minus 1 670 090-1 675 333 NA? 

Lyz2 NT_167289.2 Minus 1 707 719-1 713 452 NA? 

Lyz3 NT 167289.2 Plus 1 722 393-1 728 285 NA? 

Lyz4 NT 167289.2 Plus 1 760 703-1 766 608 NA? 

Bottlenose dolphin (Tursiops truncates) 

Lyz scaffold 114746 Plus 182 136-187 936 ENSTTRG00000013948 

Minke whale (Balaenoptera 

acutorostrata scammoni) 

Lyz NW_006733011 Minus 27 704 269-27 709 850 XM_007195043 

Killer whale (Orcinus orca) 

Lyz NW_004438568 Plus 1 319 177-1 324 490 XM 004281877 

Yangtze River dolphin (Lipotes vexillifer) 

Lyz NW. 006790307 Minus 1 455 813-1 461 420 XM 007463554 

Sperm whale (Physeter catodon) 

Lyz NW. 006716048 Minus 6 880-11 985 XM 007118874 

Dog (Canis lupus familiaris) 

Lyz 10 Plus 11 346 500-11 350 639 ENSCAFG00000000426 

Panda (Ailuropoda melanoleuca) 

Lyz GL192893.1 Plus 308 956-313 372 ENSAMEG0000001 1820 

Horse (Equus caballus) 

Lyz 6 Plus 84 276 158-84 280 173 ENSECAG00000018113 

Rhinoceros (Ceratotherium ^ simum 

simum) 

Lyz JH767750.1 Plus 25 463 742-25 467 903 ENSP00000261267_1 


a — Not annotated as a gene in Ensembl or NCBI. 
> _ Possibly missing due to incomplete gene. 
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Figure 3 Organization of lysozyme c genes in diverse Artiodactyls and relatives 
Schematic of the genomic arrangement of lysozyme c genes, and their neighbors, derived from genomic sequences in the Ensembl and NCBI databases. 
Sizes of genes, and distances between genes are not to scale. Genes shown above the lines are encoded by the plus strand, while those below are on 


the minus strand. Genomic sequences are listed in Table 3. 


2012) zebu, Bos indicus (Canavez et al, 2012); water 
buffalo, Bubalus bubalis; Tibetan antelope (chiru) (Ge et al, 
2013), Pantholops hodgsonii, and goat, Capra hircus (Dong 
et al, 2013)) are available in the NCBI database. The 
genomes of all pecoran ruminant species contained multiple 
lysozyme c genes (Table 3 and Figures 3 and S1), in accord 
with previous results (Irwin & Wilson, 1989; Irwin et al, 1989, 


2011). For most pecoran species, lysozyme c genes could 
be mapped to large genomic contigs, or chromosomes, that 
show organizations similar to that seen in the cow (Table 3 
and Figure 3). In the sheep, one gene (Lyz3) was not 
mapped to chromosome 3, but instead to an unmapped 
contig (Table 3). Since the goat genes all map to one contig 
(Table 3) it is possible that the sheep Lyz3 gene has been 
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misplaced (Figure 3), although movement to a new location 
through recombination cannot be excluded. The yak 
lysozyme c genes map to two contigs, with one containing a 
large gap that corresponds to the location where one or 
more missing lysozyme c genes might exist (Table 3 and 
Figure 3). The lysozyme c genes in both the Tibetan 
antelope and water buffalo map to two genomic contigs that 
might be adjacent in their genomes (Table 3 and Figure 3). 
Lysozyme c genes in the zebu are each on separate contigs, 
but could be arranged as seen in the cow and other pecoran 
species (Table 3 and Figure 3). 


Mosaic evolutionary histories for exons of cow lysozyme 
c genes 

To examine the evolutionary history of the cow lysozyme c 
genes, a phylogeny of the sequences was established. 
Phylogenetic trees were constructed for each exon of the 
lysozyme c genes (Figure 4) as previous analyses 
suggested that they might have experienced different 
histories (Irwin, 2004; Irwin & Wilson, 1990; Irwin et al, 1993, 
1996; Wen & Irwin, 1999). As shown in Figure 4, different 
phylogenies were identified for each exons, with similar 
trees found if different outgroup species were used or if 
phylogenies were constructed using distance or parsimony 
methods or if only synonymous substitutions were used 


(results not shown). Some consistent phylogenetic patterns 
were observed across all exons, such as the clustering of 
the Lyz2a, Lyz2b, and Lyz2c genes and Lyz3a or Lyz3c 
being closest to Lyz3b (Figure 4). In contrast, the placement 
of some genes differed greatly between exons, such as the 
placement of Lyz1 or Lyz4 (Figure 4). To test whether there 
were Statistically significant differences between the tree 
topologies estimated by each exon, we used Tree-puzzle 
(Strimmer & von Haeseler, 1996) to compare the four 
separate exon tree topologies with data for each exon. 
Despite the short lengths of some exons, at least two of the 
three alternative topologies could be excluded by all three of 
the KH, SH, and ELW statistical tests used by Tree-puzzle, 
with all three being excluded by at least one of the tests 
(Table 4). We cannot exclude the possibility that exons 2 
and 3 share an identical evolutionary history, as these trees 
were not excluded by all three of the statistical tests, but 
exons 1 and 4 have evolutionary histories that are 
incompatible with each other and with exons 2 and 3 
indicating that at least three different histories are 
represented by these four exons (Table 4). The differences 
in the topologies are unlikely to be due to convergent 
evolution acting on the lysozyme c protein sequences as the 
differences in the topologies were also seen when only 
synonymous differences were examined (results not shown). 


Table 4  Phylogenies predicted from different cow lysozyme c gene exons are significantly different 





Tree/Data Log L Difference SE KH? SH? ELW? 

Exon 1 Tree 

Exon 1 Data —703.05 0.00 1.0000 1.0000 0.9425 BEST 

Exon 2 Data —738.42 35.37 13.1770 0.0080 0.0110 0.0001 EXCLUDED? 
Exon 3 Data -722.37 19.32 9.2124 0.0270 0.1170 0.0095 

Exon 4 Data -733.42 30.37 14.3425 0.0280 0.0320 0.1540 EXCLUDED” 
Exon 2 Tree 

Exon 1 Data —844.36 113.94 22.9026 0.0000 0.0000 0.0000 EXCLUDED? 
Exon 2 Data -730.42 0.00 1.0000 1.0000 0.9438 BEST 

Exon 3 Data —746.39 15.97 8.3002 0.0400 0.3250 0.0136 

Exon 4 Data —829.74 99.32 20.0394 0.0000 0.0000 0.0000 EXCLUDED? 
Exon 3 Tree 

Exon 1 Data —391.39 59.87 13.7897 0.0000 0.0000 0.0000 EXCLUDED? 
Exon 2 Data —343.01 11.49 6.7348 0.0580 0.3190 0.0136 

Exon 3 Data —331.52 0.00 1.0000 1.0000 0.5948 BEST 

Exon 4 Data —9365.58 34.06 13.4375 0.0070 0.0160 0.0002 EXCLUDED? 
Exon 4 Tree 

Exon 1 Data —1804.06 95.06 18.8423 0.0000 0.0000 0.0000 EXCLUDEDb 
Exon 2 Data —1780.35 71.35 20.9871 0.0010 0.0020 0.0000 EXCLUDED? 
Exon 3 Data —1767.23 58.23 20.2572 0.0020 0.0070 0.0000 EXCLUDED? 
Exon 4 Data —1709.00 0.00 1.0000 1.0000 0.9756 BEST 


a — Probability of observing the tree, given the data, from the statistical one sided KH test based on pairwise SH tests (KH), the Shimodaira- 
Hasegawa test (SH), and the expected likelihood weight test (ELW) from Tree-puzzle (Strimmer & von Haesler, 1996). 


bs Probability that the data is compatible with the tree is less than 0.05. 


10 WWW.ZOOres.ac.cn 


















































































































66 Cow Lyz2a 28 Cow Lyz2a 
Cow Lyz2c 80 Cow Lyz2c 
Cow Lyz2a 89 L_ Cow Lyz2b 
Cow Lyzl 93 Cow Lyz4 
Cow Lyz4 Cow Lyz6 
100 [— Cow Lyz3b 22 Cow Lyz7 
— Cow Lyz3c i Cow Lyz5 
Cow Lyz8 Cow Lyz8 
Cow Lyz9 7 Cow Lyz10 
Cow Lyz10 Cow Lyzl 
Cow Lyzó o s po Lyz3b 
61[— Cow Lyz5 Cow Lyz3c 
IE 7027 98 Dog Lyz 
Dog Lyz E Lyz 
Panda Lyz 
002 
"002. 
5 75 [— Cow Lyz2a R 42, —Cow Lyz2b 
52 Cow Lyz2c 
064 Cow Lyz2a 
86 L- Cow 7021 
Cow Lyz8 52 Cow Lyz3a 
Cow Lyz6 šj a Cow Lyz3b 
Cow Lyz5 Cow Lyz4 
75 Cow Lyz7 Cow Lyz9 
Cow Lyz3a Cow Lyz7 
Cow Lyz3b Cow Lyz5 
—— Cow 7021 Cow Lyz6 
Cow Lyz10 55b... Cow Lyz8 
Dog Lyz Cow Lyz10 
Panda Lyz Dog Lyz 
Panda Lyz 





0.05 


Figure 4 Phylogeny of cow lysozyme c genes derive from sequences of (A) exon 1, (B) exon 2, (C) exon 3, and (D) exon 4 
Phylogenies for each of the 4 exons of the lysozyme c genes were estimated using maximum likelihood, as implemented in Mega6.06 (Tamura et al, 


2013), using the Kimura 2-paramater model with a gamma distribution, which was the best fitting model for the sequence data. Similar results were 


obtained with the neighbor-joining method or parsimony, or the use of different outgroups. Phylogenies were generated from 152, 156, 74, and 306 


aligned bases present in all sequences for exons 1, 2, 3, and 4, respectively. The presented phylogenies were bootstrapped 500 times. 


These results are in agreement with previous conclusions 
that lysozyme c genes expressed in the abomasum of 
ruminants have experienced mosaic evolution due to gene 
conversion occurring between the coding exons (Irwin, 
2004; Irwin & Wilson, 1990; Irwin et al, 1993, 1996; Wen & 
Irwin, 1999), and suggested that the 3' untranslated (exon 
4) sequences likely best reflect the evolutionary history of 
the divergent genes, as this sequence appears to have 
experienced the fewest number of concerted evolution 
events. 


Origin and evolutionary history of lysozyme c genes in 
ruminant genomes 

To better examine the evolution of the duplicated lysozyme c 
genes in ruminant species, a phylogenetic tree was 
established for the lysozyme c sequences from the diverse 
ruminants (e.g., cow, sheep, and Tibetan antelope) and their 
close relatives (e.g., pig, cetaceans, and carnivores) (Figure 
5). Exon 4 sequences were chosen to construct this 
phylogeny as they likely best reflect the divergence of the 
genes, and have experienced lower levels of concerted 
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evolution (see above). The phylogeny shown in Figure 5 
was derived by maximum likelihood, and similar phylogenies 
were generated when neighbor-joining or parsimony was 
used (results not shown). The exon 4 phylogeny shown in 
Figure 5 of the lysozyme c genes yield strong evidence for 
the orthology of 8 of 10 types of lysozyme c genes found in 
ruminants (Figure 5). Lyz3, Lyz4, Lyz5, Lyz6, Lyz7, Lyz8, 
Lyz9, and Lyz10 orthology groups each have high (88%- 
100%) bootstrap support, with the species relationships 
within each group in general accord with the accepted 
species relationships (Figure 5). This observation implies 
that these 8 genes existed in the common ancestor of 
pecoran ruminants. The phylogenetic analysis did not 
resolve Lyz1 or Lyz2 genes as monophyletic groups, but 
instead suggested some intermixing of these genes (Figure 
5). Lyz1 and Lyz2 sequences from species of tribe Bovine 
(cow, yak, zebu, and water buffalo) formed a moderately 
supported monophyletic group that had a primary 
divergence between the Lyz1 and Lyz2 sequences. The 
tribe Bovini Lyz1 and Lyz2 sequences were then grouped 
with Lyz1 sequences from the other pecoran ruminants 
(Tibetan antelope, goat and sheep), with the Lyz2 
sequences from these same species being the outgoup to 
all of the Lyz1 and Lyz2 sequences. While it is possible that 
this distribution could be explained by an ancestor having 
four genes, and pairs of genes being lost in each species, 
an alternative explanation is that the pecoran ancestor 
possessed two genes, and that a concerted evolution event 
transferred sequences from the tribe Bovini Lyz1 exon 4 
sequence to the tribe Bovini Lyz2 gene, resulting in the 
grouping of these sequences. Support for the monophyly of 
the Lyz1 and Lyz2 genes was found from phylogenies of 
exon 2 and exon 3 sequences (results not shown). These 
results suggest that the ancestor of pecoran ruminants 
possessed 10 lysozyme c genes. 

While the ancestor of modern pecoran ruminants may 
have had 10 lysozyme c genes, several extant species have 
a higher number of genes, such as cow with 14 genes and 
the zebu with 12 genes (Tables 1 and 3). The increased 
numbers of lysozyme c genes in some ruminant species 
appear to be due to lineage-specific gene duplications. The 
phylogeny presented in Figure 5 implies lineage-specific 
duplications in three genes, Lyz2, Lyz3 and Lyz7, all of 
which occurred in species (cow and zebu) of the genus Bos. 
Both cow and zebu have three Lyz2 genes (Lyz2a, Lyz2b, 
and Lyz2c) (Tables 1 and 3). Only a single Lyz2 gene was 
found in the yak, however a gap in the genome assembly 
was found at this location (Table 3), thus it is possible that 
additional Lyz2 genes exist in this genome. Better assembly 
of the Bos genome sequences are needed to determine 
whether the triplicated genes have a single origin, or 
represent parallel duplication, a conclusion that does have 
some support from the phylogenetic analysis (Figure 5). 
Duplicated Lyz3 genes were only found in the cow, although 
the lack of this gene in the yak, potentially due to a gap in 
the assembly, and the poor assembly of the zebu genome 
do not rule out the possibility that multiple Lyz3 genes exist 
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in these species (Tables 1 and 3). The duplications of the 
Lyz2 and Lyz3 genes in Bos likely represent products of 
segmental duplications (Liu et al, 2009; Seo et al, 2013). It is 
possible that segmental duplications may also exist in other 
pecoran ruminant species, but were collapsed as single 
genes during the genome sequence assembly process, and 
thus the increased numbers seen in the genus Bos simple 
reflect the better cow genome assembly. 

The distribution of the numbers of lysozyme c genes in the 
genomes of ruminant species and their close relatives is 
consistent with an amplification of the lysozyme c gene on 
the lineage leading to true ruminants as previously proposed 
(Irwin & Wilson, 1989; Irwin et al, 1989, 1992, 2011; Yu & 
Irwin, 1996). In contrast, our current phylogenetic analysis of 
the 3' untranslated regions of lysozyme c genes suggests 
that the amplification of these genes was initiated very early 
in the artiodactyl lineage, before the divergence of the rumi- 
nants and tylopod (e.g., alpaca) lineages, and implying that 
the pig and cetaceans have lost genes, however these early 
divergences are very poorly supported (Figure 5). Indeed, 
phylogenetic analysis of exon 1, exon 2, or exon 3 
sequences by themselves yielded differing conclusions 
concerning these earliest duplications, although again, none 
of these analysis yielded strong conclusions (results not 
shown). Analysis of larger amounts of genomic sequences 
(e.g., intronic and flanking sequence) potentially could 
resolve the order of the earliest divergences of the 
paralogous lysozyme c genes and cetartiodactyl species. 
While the alpaca has multiple lysozyme c genes (Table 3), 
our phylogenetic analysis suggests that they originated 
through a parallel series of lineage-specific independent 
duplications. 


Rates of evolution in ruminant lysozyme c genes 

Duplication of the lysozyme c gene on the ruminant lineage 
has allowed the specialization of gene expression in distinct 
tissues, such as different chambers of the stomach, and 
thus evolution of novel gene function (Callewaert & Michiels, 
2010; Jiang et al, 2014; Irwin et al, 1992; Irwin, 1995, 2004; 
Prager & Jollés, 1996). Changes in the function of lysozyme 
c likely leads to changes in the evolutionary constraints 
acting upon these genes. To examine this issue we 
calculated the divergence at nonsynonymous and 
synonymous sites among lysozyme c genes, with the results 
from three divergent representatives of pecoran ruminants 
(cow, goat, and Tibetan antelope), and between genes in 
these three species and the single copy lysozyme c gene 
sequences found in pig and horse shown in Table 5 (similar 
results were seen with the other pecoran ruminant species). 
The relative rates of nonsynonymous to synonymous 
substitutions (dn/ds) varied between genes when compared 
among ruminants, from low values for the Lyz5 and Lyz6 
genes, which imply that they are strongly constrained, to 
high values for the Lyz3 and Lyz9 genes, suggesting that 
there is little constraint on their protein sequences (Table 5). 
The cow Lyz3 and Lyz9 genes fail to predict intact open 
reading frames, suggesting that they are pseudogenes 
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Figure 5 Phylogeny of ruminant lysozyme c genes derived from exon 4 sequences of predicted lysozyme c genes 

The phylogeny of the lysozyme c genes was estimated from aligned exon 4 sequences (192 aligned bases in all sequences) using maximum likelihood, 
as implemented in Mega6.06 (Tamura et al, 2013), using the Kimura 2-paramater model with a gamma distribution, which was the best fitting model for 
the sequence data. Similar results were obtained with the neighbor-joining method or parsimony. The phylogeny was bootstrapped 500 times. Outgroups 
used to root the phylogeny are shown at the bottom. The ten types of lysozyme c genes are indicated on the right, with the bootstrap values that support 
8 of these clades (all except the Lyz1 and Lyz2 clade) shown in bold. 
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Table 5 Rates of evolution of ruminant lysozyme c genes 




















Lyz1 Lyz2 Lyz3 Lyz4? Lyz5 

dn ds dn/ds dn ds dn/ds dn ds dn/ds dn ds dn/ds dn ds dn/ds 
Cow-Goat 0.026 0.043 0.611 0.052 0.081 0.640 0.073 0.054 1.350 0.031 0.048 0.647 0.022 0.114 0.197 
Cow-Tibetan antelope 0.041 0.033 1.234 0.045 0.058 0.767 0.043 0.032 1.340 NA NA NA NA NA NA 
Goat-Tibetan antelope 0.026 0.043 0.610 0.007 0.066 0.105 0.028 0.021 1.314 NA NA NA NA NA NA 
Pecora average 0.818 0.504 1.334 0.647 0.197 
Pig-Cow 0.184 0.327 0.564 0.246 0.342 0.720 0.268 0.268 0.998 0.230 0.264 0.872 0.224 0.354 0.632 
Pig-Goat 0.192 0.317 0.604 0.231 0.318 0.726 0.271 0.259 1.048 NA NA NA NA NA NA 
Pig-Tibetan antelope 0.193 0.272 0.710 0.233 0.307 0.758 0.254 0.287 0.886 NA NA NA NA NA NA 
Pig average 0.626 0.735 0.977 0.872 0.632 
Horse-Cow 0.099 0.438 0.226 0.180 0.467 0.386 0.129 0.459 0.281 0.193 0.345 0.560 0.197 0.394 0.501 
Horse-Goat 0.083 0.451 0.184 0.171 0.441 0.388 0.144 0.379 0.380 NA NA NA NA NA NA 
Horse-Tibetan antelope 0.101 0.372 0.272 0.171 0.482 0.356 0.127 0.404 0.315 NA NA NA NA NA NA 
Horse average 0.227 0.377 0.325 0.560 0.501 

Lyz6 Lyz7 Lyz8 Lyz9 Lyz10 

dn ds dn/ds dn ds dn/ds dn ds dn/ds dn ds dn/ds dn ds dn/ds 
Cow-goat 0.025 0.128 0.195 0.043 0.098 0.441 0.058 0.067 0.867 0.052 0.047 1.097 0.033 0.058 0.573 
Cow-Tibetan antelope 0.019 0.073 0.258 0.049 0.078 0.627 0.054 0.083 0.652 0.038 0.031 1.215 0.046 0.072 0.635 
Goat-Tibetan antelope 0.025 0.065 0.385 0.049 0.071 0.692 0.023 0.028 0.825 0.029 0.014 1.996 0.012 0.026 0.446 
Pecora average 0.279 0.587 0.781 1.436 0.551 
Pig-Cow 0.215 0.411 0.523 0.229 0.382 0.599 0.256 0.389 0.658 0.247 0.352 0.701 0.191 0.270 0.708 
Pig-Goat 0.238 0.427 0.557 0.278 0.445 0.625 0.265 0.316 0.840 0.272 0.393 0.691 0.187 0.248 0.753 
Pig-Tibetan antelope 0.218 0.351 0.621 0.261 0.386 0.677 0.260 0.330 0.788 0.253 0.392 0.645 0.190 0.271 0.700 
Pig average 0.567 0.634 0.762 0.679 0.721 
Horse-Cow 0.185 0.432 0.429 0.193 0.419 0.461 0.221 0.449 0.491 0.129 0.368 0.349 0.113 0.351 0.323 
Horse-Goat 0.204 0.411 0.497 0.227 0.457 0.498 0.241 0.416 0.579 0.169 0.406 0.415 0.101 0.366 0.277 
Horse-Tibetan antelope 0.189 0.366 0.516 0.220 0.443 0.498 0.230 0.434 0.531 0.152 0.400 0.380 0.094 0.363 0.259 
Horse average 0.481 0.485 0.534 0.381 0.286 


? — Sheep Lyz4 used to replace the incomplete goat Lyz4 for comparisons. 


(Table 1 and Figure 2) and thus should have no evolutionary 
constraints on their protein sequences. 

The cow genome contains three Lyz3-like genes, with only 
one being a full-length gene sequence (Lyz3b), and a single 
copy of this gene was found in most of the other pecoran 
ruminant species (Tables 1 and 3). Cow Lyz3b gene was 
previously identified as the cow lysozyme c 
pseudogene NS4 (Irwin, 1995, 2004). While the Lyz3 gene 
sequences from tribe Bovini (cow, zebu and water buffalo) 
all share a frame shift mutation in exon 3 (amino acid 
residue 100 in Figure 2), which would prevent translation of 
a functional product, and additional mutations that 
potentially disrupt functions found in some sequences at 
other locations, the sequences from sheep, goat, and 
Tibetan antelope all predict a full-length open reading frame 
(Figure S1). This observation might suggest that the Lyz3 
gene became a pseudogene, due to a frame-shifting 
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mutation, on the lineage leading to tribe Bovini, after 
divergence from the other pecoran ruminant lineages. 
However, a high rate of nonsynonymous substitutions is 
also observed between the goat and Tibetan antelope Lyz3 
gene sequences (Table 5) and between the sheep and 
both the goat and Tibetan antelope sequences (results not 
shown) suggesting that few evolutionary constraints were 
acting on this sequence and that this gene may have been 
non-functional in the common ancestor of all pecoran 
ruminants. It is possible that a mutation that prevented 
expression, or an amino acid substation that that preve- 
nted function, rather than a mutation that prevents 
translation of an intact product, was the initial mutation that 
created this pseudogene. 

The second gene with a very high dn/ds ratio is the Lyz9 
gene, which is composed of only 2 exons, exons 1 and 4, in 
the cow, due to the loss of exons 2 and 3 (Table 1 and 


Figure 2). Orthologs of the Lyz9 genes in other pecoran 
ruminant species also have similar gene structures (Tables 1 
and 3 and Figure 51), suggesting that this structure exited in 
the Lyz9 gene in the ancestor of all pecoran ruminants. The 
loss of exon 2 and 3 sequences from Lyz9 prevents the 
translation of a functional lysozyme, thus it can be 
concluded that this pseudogene originated before the 
radiation of the pecoran ruminants. Consistent with this 
conclusion, a high rate of divergence at nonsynonymous 
Sites is observed in the Lyz9 gene sequence among all 
pecoran ruminant species (Table 5 and results not shown). 

The Lyz3 and Lyz9 genes account for 4 of the five 
predicted lysozyme c pseudogenes in the cow genome 
(Table 1 and Figure 2). In addition to a pair of inframe stop 
codons (located between amino acid residues 24 and 25, 
and residue 26, Figure 2), the initiation codon for the cow 
Lyz8 gene is valine rather than methionine (amino acid -18 
in Figure 2 and Figure S1). Orthologs of the Lyz8 gene from 
members of the tribe Bovini (yak, zebu, and water buffalo) 
share the inframe stop codons, as well other mutations such 
as a 9 base deletion in exon 2 (3 codons - residues 66-68 in 
Figure 2), while Lyz8 sequences from other pecoran species 
(sheep, goat and Tibetan antelope) do not possess any 
obvious harmful amino acid substitution, other than the 
valine substitution at the initiation codon (Figure S1). In 
contrast to the Lyz3 and Lyz9 pseudogenes, a much lower 
dn/ds ratio was observed in the pairwise comparisons 
among cow, goat, and Tibetan antelope (Table 5), which 
would be consistent with functional constraints acting on 
some, but not necessarily all, of the Lyz8 protein sequences. 
These observations appear to suggest, that despite the 
replacement of the initiator methionine with valine, the Lyz8 
protein sequences in the sheep, goat and Tibetan antelope 
is functional, while a mutation occurred on the lineage 
leading to tribe Bovini to producing the Lyz8 pseudogene. 
How a functional protein can be translated from the Lyz8 
gene, or evolutionary constraints that mirror protein function, 
is unclear. A downstream ATG, at codon 85 of the mature 
protein sequence (Figure 2), would be predicted to yield a 
protein of only 45 amino acid residues, far shorter than the 
typical 145 amino acid long protein lysozyme c precursor, 
with most of the sequence not being translated and thus not 
under evolutionary constraint for protein function. 


Episodic evolution of ruminant lysozyme c genes 

The cow lysozyme c genes displaying the lowest dn/ds 
ratios among ruminant species, and thus implying the 
strongest evolutionary constraints, are the Lyz5 and Lyz6 
genes (Table 5), which are expressed in the abomasum 
(Table 1). Lysozyme c genes expressed predominantly in 
non-stomach tissues (Irwin, 2004), such as Lyz1, in milk, 
Lyz2, in the trachea, and Lyz10, in macrophages, have 
intermediate dn/ds ratios, but ratios that lower than those 
seen for the Lyz3, Lyz8, and Lyz9 pseudogenes (Table 5). 
However, when the dn/ds ratios are calculated between the 
ruminant genes (cow, goat, and Tibetan antelope) and an 
outgroup sequence (pig or horse), the stomach expressed 


Lyz5 and Lyz6 genes are seen to have dn/ds ratios that are 
either similar (when pig is the outgroup) or higher (horse 
being the outgroup) than those seen for the non-stomach 
(Lyz1, Lyz2, and Lyz10) genes (Table 5). To obtain this 
pattern of results, these observations suggest that the dn/ds 
ratio on the common ancestral lineage leading to the 
ruminants, after divergence from pig or horse, but before 
radiation of the pecoran ruminants, was higher for the 
lysozyme c genes expressed in the abomasum (Lyz5 and 
Lyz6) than for those expressed in non-stomach tissues 
(Lyz1, Lyz2, and Lyz10). This suggests that the rates of 
evolution of lysozyme c genes expressed in the abomasum 
display an episodic pattern, with more rapid evolution on the 
early ruminant lineage, and a slower rate within the pecoran 
ruminants. These results are consistent with previous 
findings of accelerated evolution of lysozyme c protein 
sequences obtained from the abomasum of ruminant 
species (Jollés et al, 1989; Irwin & Wilson, 1990; Irwin et al, 
1992, 1993). 


CONCLUSIONS 


Genome sequences have advanced our understanding of 
the evolution of the lysozyme c gene family in ruminant 
species. Genomic sequences from seven divergent pecoran 
ruminant species allowed us to demonstrate that the 
genome of the pecoran ruminant common ancestor 
possessed at least 10 lysozyme c genes, and that these 
genes have largely been retained by extant ruminant 
species. More recent gene duplication, likely via segmental 
duplications (Liu et al, 2009; Seo et al, 2013), have resulted 
in increases in the number of lysozyme c genes on some 
lineages, with 14 genes found in the cow, but we can not 
exclude the possibility that some duplications may have 
been missed during assembly of some genomes. Lysozyme 
c genes have not evolved in a simple divergent manner, but 
rather by concerted evolution acting independently on each 
exon, yielding differing phylogenetic relationships for the ten 
types of lysozyme c genes. Some lysozyme c genes have 
become pseudogenes, either due to mutations in their 
coding sequence (e.g., Lyz3 and Lyz8) or by deletion of 
exon sequences (e.g., Lyz9). Some pseudogenes may have 
been generated by incomplete duplication of genes, such as 
Lyz3a and Lyz3c in the cow. Despite being presumably non- 
functional, at least two pseudogenes that exited in the 
ancestral pecoran ruminant (Lyz3 and Lyz9) have been 
retained in diverse ruminant species. A third lysozyme c 
gene has its initiation codon mutated to valine (from 
methionine), yet shows evidence that its coding sequence is 
evolutionary constrained on some ruminant lineages. This 
suggests that some lysozyme c pseudogenes may retain 
biological functions, however, how protein function in this 
sequence is maintained is unclear. Changes in the rates of 
nonsynonymous substitutions suggest that changes have 
occurred in the functional constraints acting on lysozyme c 
protein sequences, and these changes have occurred in an 
episodic fashion. 
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ABSTRACT 


The coexistence of ecologically similar species 
sharing sympatric areas is a central issue of 
community ecology. Niche differentiation is required 
at least in one dimension to avoid competitive 
exclusion. From 2012-2014, by adopting the 
methods of mist-nets and point counts to evaluate 
spatial niche partitioning and morphological 
differentiations, we explored the coexistence 
mechanisms of seven sympatric fulvettas in Ailao 
Mountains, Ejia town, Yunnan Province, China. The 
microhabitats of these seven fulvettas were 
significantly different in elevation, roost site height 
and vegetation coverage, indicating a spatial niche 
segregation in different levels. Approximately, 
90.3096 of the samples were correctly classified by 
linear discriminant analysis (LDA) with correct rates 
at 91.20%-100%, except the White-browed fulvetta 
(Alcippe vinipectus) (65.4%) and the Streak-throated 
fulvetta (A. cinereiceps) (74.696). The seven fulvettas 
were classified into four guilds based on their 
specific morphological characters, suggesting that 
the species in each guild use their unique feeding 
ways to realize resource partitioning in the 
overlapped areas. These finding indicate that 
through multi-dimensional spatial niche segregation 
and divergence in resource utilizing, the inter- 
specific competition among these seven fulvettas is 
minimized, whereas, coexistence is promoted. 


Keywords:Fulvettas;Coexistence;Niche segregation; 
Ailao Mountains 


INTRODUCTION 
The competition and coexistence of ecologically similar species 


sharing sympatric areas is one of the hot topics in community 
ecology (Macarthur, 1958). According to the competitive 
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exclusion principle, to reduce inter-specific competition, niche 
segregation is required at least in one dimension among 
sympatric congeners (Bagchi et al, 2003; Denoél et al, 2004; 
Munday et al, 2001; Schoener, 1974; Tschapka, 2004). The 
differentiation in resource use can be taken as a standard to 
judge niche partitioning (Fox, 2004; Svenning, 1999). Lots of 
studies have shown that niche partitioning was principally along 
three dimensions: time, space and diet (Chesson, 2000; Davies 
et al, 2007; Martinez-Freiria et al, 2010; Schoener, 1974). 
Temporal partitioning includes daily (Di Bitetti et al, 2009; Lara 
et al, 2011; Lucherini et al, 2009) or seasonal differences 
(Martínez-Freiría et al, 2010; Schuett et al, 2005) in animals' 
activity patterns. Spatial partitioning includes differences in 
habitat selection (Chiang et al, 2012; Quillfeldt et al, 2013) or in 
microhabitat utilization (Langkilde & Shine, 2004; Vidus-Rosin 
et al, 2012). Trophic partitioning includes differences in prey 
size (Kaifu et al, 2013) or prey type (Ward-Campbell et al, 2005). 
Niche segregation may occur along either some combination of 
these dimensions or just upon one of them (Loveridge & 
Macdonald, 2003). Furthermore, niche segregation could be 
related to species’ morphology (Gustafsson, 1988) because 
different morphological characteristics may result in different 
behaviors (Miles et al, 1987) and different behaviors allow 
animals to partition the limited resources in different ways, 
hence reduce the inter-specific competitive interactions 
(Guillemain et al, 2002). Therefore, it is necessary to evaluate 
the influences of various factors in different ways as well as at 
different scales when exploring the mechanisms of animals’ 
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coexistence (Li et al, 2013). 

Closely related species usually have similar ecological 
characteristics and occupy similar ecological niches. When 
they contact, the inter-specific competition may arouse (Wiens 
& Graham, 2005). Studies on the sympatric congeners of 
birds in China were focused on species such as pheasants 
(Cui et al, 2008; Li et al, 2006), herons (Wen et al, 1998; Ye et 
al, 2006; Zhu et al, 1998), tits (Gao & Yang, 1991; Liu et al, 
1989; Yang et al, 2012), woodpeckers (Gao et al, 1997), larks 
(Zhao & Zhang, 2004) and prinias (Zhou & Fang, 2000), but 
never on fulvettas. Alcippe is a group of Timaliidae includes 
18 species worldwide (Zheng, 2002) with 15 of them 
distributing in China (Zheng, 2011) and 12 of them in Yunnan 
Province, in specific (Yang & Yang, 2004). The range of body 
length in fulvettas is 11-15 cm. Different species of fulvettas 
have similar morphology, plumage color and stripes and 
sexual differences are difficult to distinguish morphologically. 
They inhabit the undergrowth in broadleaf forests, mixed 
coniferous broad-leaved forests, bamboo groves, scrubs and 
brambles and feed on animal-based food, such as insects, 
caterpillars and mollusks, but sometimes on plant food (Yang 
& Yang, 2004). The Golden-breasted Fulvetta (Alcippe 
chrysotis), Rufous-winged Fulvetta (A. castaneceps), White- 
browed Fulvetta (A. vinipectus), Spectacled Fulvetta (A. 
ruficapilla), Streak-throated Fulvetta (A. cinereiceps), Rusty- 
capped Fulvetta (A. dubia) and Grey-cheeked Fulvetta (A. 
morrisonia) are the seven species of fulvettas coexisting in 
Ailao Mountains, Ejia town, Yunnan Province, China, which are 
excellent subjects in the study on coexistence mechanisms of 
sympatric congeners. In this study, we aim to: (1) understand 
the coexistence mechanisms of these fulvettas through 
exploring their spatial niche partitioning and morphological 
differentiation; (2) test the niche theory on animal communities 
in the subtropical mountain forest; and (3) discuss the necessity 
to evaluate the influences of various factors in different ways as 
well as at different scales when exploring the coexistence 
mechanisms of sympatric congeners. 


MATERIALS AND METHODS 


Study area 

The study site consisted of two adjacent areas, the Ailaoshan 
National Nature Reserve and Konglong River Nature Reserve, 
was located in the steep (30°- 40°, or even 60°) eastern slope 
(800-2 800 m a.s.l.) of north central of Ailao Mountains, Ejia 
Town, Yunnan Province (Chen & Ye, 1988; Liu et al, 1988; Wu 
& Yang, 2008). Because of the west monsoon climate, annual 
temperature differences are small, whereas, daily temperature 
differences are large. The large elevation gradients also lead 
to the vertical climatic spectrum. The radiation, rainfall, tempe- 
rature and forest types are all featured with obvious vertical 
variations (Wang et al, 1988). Four major forest types from 
low to high elevation are observed: (1) savanna shrub and 
grass; (2) dry evergreen broadleaved forest; (3) semi-moist 
evergreen broadleaved forest and Burma pine (Pinus 
yunnanensis) forest; and (4) moist evergreen broadleaved 
forest. However, most of the original vegetation in lower 


elevations have been destroyed or replaced by cultivated land 
and villages due to long-term human disturbance (Liu et al, 
1988). 


Birds sampling 

Fieldwork was conducted during breeding seasons and 
winters of 2012-2014. Mist-nets work best in field survey 
because fulvettas inhabit the undergrowth (Bibby et al, 1998). 
Point counts are effective in species survey and in montane 
forest bird community survey (Wu et al, 2010). Mist-nets 
combined with point counts will offset the weakness of each 
other, and provide a more accurate survey result (Pagen et al, 
2002; Rappole et al, 1998). Birds inhibit in lower elevations 
start to reproduce earlier due to the higher temperature there. 
So, in breeding seasons, the fieldwork was carried out from 
valley to montane crest, whereas, in winters, the order was 
inversed. 

We divided the study area into 10 units along elevation 
gradients with each unit containing a 200 m elevation 
differences. Ten mist-nets of the same specification (12 mx2.5 
m, 36 mm mesh) were set up symmetrically in each unit. Due 
to the rugged terrain and limited access, it was difficult to set 
up each mist-net with 20 m elevation differences exactly. 
When different habitat types occurred in one unit, mist-nets 
were set up in each different habitat proportionally. The 
location, elevation and working time were recorded by GPS 
(NAVA 100) and the habitat type as well as environment 
information were also recorded meanwhile. We chose a 
rectangular patch (with the diagonal of 20 mx20 m) around 
the mist-net to estimate vegetation coverage. Observers 
walked along the diagonal of the rectangular patch to record 
the vegetation (tree, shrub and herb) of every 1 m. At each 
point of 1 m, if vegetations were found, then it was recorded 
as 1, while, if not, it was 0. Then the vegetation coverage of 
each type (tree, shrub and herb) was estimated according to 
the percentages of the vegetations recorded at all points. 
Mist-nets were kept open for 3 days and were remained 
closed during raining. We checked each mist-net hourly 
during the day and marked each captured individual with 
metal ring. Time, species and numbers of birds were recorded. 
Body weight (BW) was measured by an electronic balance 
(DIAMOND,  precisionz0.1 g) and other morphological 
characteristics such as body length (BL), wing length (WL), 
tail length (TL), tarsus-metatarsus length (TML), claw length 
(CLL), finger length (FL) and culmen length (CUL) were 
recorded by a vernier caliper (precision=0.05 mm) according 
to Zheng (1995). All the juveniles or certain species with too 
small sampling size were excluded from measuring. 

Unlimited radius point counts were conducted in our study 
and the points were established on either pre-existing trails or 
low traffic volume roads along the elevation gradient in each 
unit. The points were systematically 200 m apart in a three- 
dimensional space. Surveys were conducted during the peak 
period of birds' activities from sunrise to 4 h after sunrise. We 
used GPS receiver (NAVA 100) to record the location of each 
point. There were 160 points in total in our study area and 
each count lasted for 10 min. During this period, we used 
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binoculars (Eagle 8x40) to observe birds and the indentified 
ones were recorded (Table 1). 


Table1 Information recorded by point counts 


Items Details 





Horizontal distance from the observer to the detected 


Distance bird or the average horizontal distance to a group of 
birds 
Position of a bird was initially sighted (A: trunk of a tree; 
" B: branches of a tree; C: substratum of crown canopy; 
Roost site ; 
D: superstratum of crown canopy; E: inner part of a 
shrub; F: outer part of a shrub; G: ground) 
,. Vertical distance from ground to the position of a bird 
Roost site ARH j iti 
heioht was initially sighted or to the central position of a group 
7 of birds 
Time Time when the bird was detected 
Animal 
위 Bird species and numbers of each species 
subject 


Data analysis 

Statistical analysis was conducted using R software (R Core 
Team, 2013). Shapiro and Levene statistics were used to 
test for normality and homogeneity. Function “pair-wise-t- 
tes" was used to conduct multiple comparisons for 
unbalanced design when variance analysis was afforded 
(Crawley, 2012). Kruskal-Wallis rank sum test was used 
when data was not in normality or homogeneity (Crawley, 
2012) and function "gao cs" in the "nparcomp" package 
(Gao et al, 2008) was applied to conduct unbalanced 
multiple comparisons. The tested variables were 
morphological characteristics and habitat variables. All the 
data was displayed in mean+SD. Statistical tests were two- 
tailed-tests and the confidence level was 95%. We used 
model selection based on generalized linear models (GLM) 
with Poisson error to find major factors affecting birds’ 
habitat selection. Morphological characteristics were applied 
in linear discriminant analysis (LDA) to estimate the 
similarity among different species and to predict resource 
partition. Morphological data was standardized by dividing 
the cubic root of body weight to avoid the influence of body 
size (Amadon, 1943). Culmen length was excluded from 
standardization because it was primarily related to the size 
of food (Hespenheide, 1973). Logarithmic transformation 
was conducted before the relative length data (L=L/BW1/ 3, 
L: the relative length of morphological characteristics) was 
applied to LDA (Z’=logioZ, Z: the relative length data) 
(Atchley et al, 1976). Data collected by mist-nets and point 
counts were analyzed independently to avoid influences of 
different methods (Barlow et al, 2007). 


RESULTS 


Habitat variables 
Elevational distribution All of the seven fulvettas in Ailao 
Mountains, Ejia Town, Yunnan Province, were documented in 
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this study and their elevation distributions were demonstrated in 
Figure 1. The Grey-cheeked Fulvetta occupied the lowest 
elevation, while the  Golden-breasted Fulvetta preferred 
relatively higher elevation. Numbers documented by mist-nets 
and point counts was 413 and 417, respectively (Table 2). Signi- 
ficant differences in elevation were detected both by mist-nets 
(y7=253.6158, P<0.01) and point counts (y7= 269.1986, P«0.01). 
Further analysis of data from mist-nets showed that there were 
no significant differences among White-browed Fulvettas, Streak- 
throated Fulvettas and Golden-breasted Fulvettas and similar 
results were also occurred between Rusty-capped Fulvettas and 
Spectacled Fulvettas. Significant differences were detected 
among the other pairs. Rufous-winged Fulvettas were excluded 
from analyses because of the small sample size (only two 
records). However, non-significant differences were only found 
between Streak-throated Fulvettas and Rufous-winged Fulvettas 
in all of the 21 combinations of point counts (Table 3). 
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Figure 1 Scatter diagram of fulvetta species’ distribution along 
elevation gradients 

1: Grey-cheeked Fulvetta (A.morrisonia); 2: Rusty-capped Fulvetta 
(A.dubia); 3: Spectacled Fulvetta (A.ruficapilla); 4: Streak-throated Fulvetta 
(A.cinereiceps); 5: Golden-breasted Fulvetta (A.chrysotis); 6: White-browed 
Fulvetta (A.vinipectus); 7: Rufous-winged Fulvetta (A.castaneceps). 


Roost site choice and roost site height Significant 
differences (x^-248.811, P«0.01) were found in the roost site 
height of the 417 fulvettas documented by point counts. Only 
three combinations consisted of White-browed Fulvettas, 
Spectacled Fulvettas and Streak-throated Fulvettas did not 
show significant differences in all of 21 combinations of the 
seven fulvettas in nonparametric multiple tests (Table 4). 
Rufous-winged Fulvettas occupied the highest roost site while 
Rusty-capped Fulvettas held the lowest one (Table 4). 
Although five species mainly chose the inner part of a shrub 
to roost, subtle differences were still detected (Figure 2). 
Golden-breasted Fulvettas preferred the substratum of crown 
canopy and Rufous-winged Fulvettas preferred branches of 
trees. White-browed Fulvettas and Streak-throated Fulvettas, 
which did not show significant differences in distribution 
elevations or roost site height, made similar choices over 
roost sites located in the inner part of shrubs, and the 
percentage was more than 80%. 


Table 2 Numbers of each fulvetta species recorded by mist-nets or point counts 


Numbers of each species recorded by mist- 


Numbers of each species recorded by point 





Species nets (n) counts (n) 
White-browed Fulvetta (A. vinipectus) 31 22 
Rufous-winged Fulvetta (A.castaneceps) 2 23 
Golden-breasted Fulvetta (A.chrysotis) 32 46 
Streak-throated Fulvetta (A.cinereiceps) 74 92 
Spectacled Fulvetta (A.ruficapilla) 85 50 
Rusty-capped Fulvetta (A.dubia) 41 78 
Grey-cheeked Fulvetta (A.morrisonia) 148 106 

Total 413 417 


Table 3 Comparisons of the distribution elevations of the fulvetta species recorded by mist-nets or point counts (m) 


Species 


Elevations of the species recorded by mist-nets Elevations of the species recorded by point counts 





Rufous-winged Fulvetta (A.castaneceps) 

2 633.742 + 108.005? 
2 600.438 + 163.993° 
2 554.919 +223.930° 
2 190.459 + 211.938" 
2 090.634 + 332.147° 
1 838.588 + 334.220* 


White-browed Fulvetta (A. vinipectus) 
Golden-breasted Fulvetta (A.chrysotis) 
Streak-throated Fulvetta (A.cinereiceps) 
Spectacled Fulvetta (A.ruficapilla) 
Rusty-capped Fulvetta (A.dubia) 





Grey-cheeked Fulvetta (A.morrisonia) 


Same superscripts indicate non-significant differences. 


2 552.130+66.213° 
2 516.455+65.236° 
2 576.478 - 31.804? 
2 559.185-- 140.016” 
2 159.380 +234.746° 
2 023.295 +361.512° 
1 728.330+260.232! 





Table 4 Sampling size and roost site height of the seven fulvetta species 





Species Numbers (n) Roost site height (m) 
Rufous-winged Fulvetta (A.castaneceps) 23 11.087+2.109° 
White-browed Fulvetta (A. vinipectus) 22 3.38616.353* 
Golden-breasted Fulvetta (A.chrysotis) 46 5.978+3.480° 
Streak-throated Fulvetta (A.cinereiceps) 92 1.328+1.345° 
Spectacled Fulvetta (A.ruficapilla) 50 1.300+0.909° 
Rusty-capped Fulvetta (A.dubia) 78 0.559+0.408° 
Grey-cheeked Fulvetta (A.morrisonia) 106 3.774+2.737° 


Same superscripts indicate non-significant differences. 


Vegetation coverage and mountain slope inclination 
Significant differences were found in tree coverage (x°=34.9663, 
P«0.01), shrub coverage (x’=49.2832, P«0.01), herb coverage 
(y’=50.828, P«0.01) as well as in the shrub and herb coverage 
between White-browed Fulvettas and Streak-throated Fulvettas 
(Table 5). Because no niche segregation was found between 
these two Fulvetta species in earlier analysis, these differences in 
vegetation coverage indicate the specific choices of these two 
species over microhabitat utilization and their relaxed inter- 
specific competition. No significant differences were found in the 
mountain slope inclination (X= 10.1526, P=0.1184). Rufous- 
winged Fulvetta was excluded from analyses due to the small 
sampling size. 


Key factors affect habitat selection 

According to earlier analyses of habitat variables, we chose 
the elevation, tree coverage, shrub coverage and herb 
coverage as habitat selection factors and used function 
“aictab” in the “AlCcmodavg” package (Mazerolle, 2013) to 
determine the key factors affecting habitat selection. The 
results showed that the first 4 models’ cumulative weight of 
QAICc (Cum.Wt) reached 100%. The first two models’ delta 
QAICc were<2 and their cumulative weight reached 81%. 
The first one had the smallest QAICc and its weight reached 
57% (Table 6). Because the weight of the first model was 
more than twice of the second one and it was succinct, the 
first model was considered as the optimal model. Because the 
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* Groud = Substratum of crown canopy 
* Outer part of shrub = Branch of trees 

* [nner part of shrub = Trunk of trees 

* Superstratum of crown canopy 


Figure 2 Comparison of roost site choices of species 

1: White-browed Fulvetta (A.vinipectus); 2: Streak-throated Fulvetta 
(A.cinereiceps); 3: Golden-breasted Fulvetta (A.chrysotis); 4: Rufous- 
winged Fulvetta (A.castaneceps); 5: Spectacled Fulvetta (A.ruficapilla); 6: 
Rusty-capped Fulvetta (A.dubia); 7: Grey-cheeked Fulvetta (A.morrisonia). 


first model was composed of elevation, herb coverage and 
shrub coverage, these factors were taken as the key factors 


affecting fulvettas' habitat selection. 

We conducted model average of these four factors 
(elevation, tree coverage, shrub coverage and herb coverage) 
to learn effect trends. The results (Table 7) also showed that 
tree coverage did not significantly affect fulvettas' habitat 
selection because 0 occurred within its 95% confidence 
interval. The estimators of elevation and shrub were > O0, 
indicating their positive effects on fulvettas' habitat selection. 
Herb coverage had negative effects on fulvettas' habitat 
selection because its estimator<0. Therefore, we concluded 
that fulvettas prefer the habitats with relatively high elevation 
and thick shrubs. 


Morphological characteristics 

Significant differences were detected in body weight 
(y7=292.4995, P<0.01), body length (y*=273.2875, P<0.01), 
wing length (x?=270.8071, P«0.01), tail length (y°=252.251, 
p<0.01), tarsus-metatarsus length (x7=183.7224, P«0.01), 
claw length (x?=81.1642, P<0.01), finger length (y7=169.0881, 
P«0.01) and culmen length (y*=264.6239, P<0.01). The 
Rufous-winged Fulvetta had the smallest tail length and 
tarsus-metatarsus length, while, the Golden-breasted Fulvetta 
was the smallest in the other morphological characteristics. 
Data of the Rusty-capped Fulvetta and the Grey-cheeked 
Fulvetta were relatively large (Table 8). 


Table 5 Comparisons of vegetation coverage of the habitats of seven fulvetta species 


Species 

i (%) 
White-browed Fulvetta (A. vinipectus) 60.461 +30.2047"° 
Golden-breasted Fulvetta (A.chrysotis) 64.8534 14.120°° 


Streak-throated Fulvetta (A.cinereiceps)  67.850+24.946 


Spectacled Fulvetta (A.ruficapilla) 53.357 + 26.087° 


Rusty-capped Fulvetta (A.dubia) 63.019+ 25.633" 





Grey-cheeked Fulvetta (A.morrisonia) 71.7193 24.159? 


Same superscripts indicate non-significant differences. 


Percentage of tree coverage 


Percentage of shrub coverage Percentage of herb coverage 


(%) (%) 
40.268 +37.980° 78.148+39.234° 
61.050+17.664° 84.141 417.565? 
72.761 +32.170° 52.393+41.027° 
80.808 + 26.292? 38.211+43.193° 


72.461 +30.159? 40.637 + 43.533° 








77.666 + 24.346? 34.364 +42.520° 


Table 6 Results of model selection according to quasi-likelihood of the second-order Akaike’s information criterion 





Candidate models QAICc Delta QAICc QAICcWt Cum.Wt Quasi.LL 
1/3/4 295.87 0.00 0.57 0.57 —142.83 
1/2/3/4 297.62 1.75 0.24 0.81 —142.67 
1/3 299.14 3.27 0.11 0.92 —145.5 
1/4 299.77 3.9 0.08 1 —145.82 
1 312.52 16.65 0 1 —153.22 
3/4 322.2 26.33 0 1 -157.03 
3 325.07 29.2 0 1 —159.49 
4 339.48 43.61 0 1 —166.7 





1: Elevation; 2: Percentage of tree coverage; 3: Percentage of shrub coverage; 4: Percentage of herb coverage. 
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Table 7 Model-averaged parameter estimations of habitat factors 






































Factors Estimator SE Lower 95% Cl Upper 95% Cl 

Elevation 0.0012 2e-04 76-04 0.0016 

Percentage of tree coverage —0.0021 0.0037 —0.0094 0.0051 

Percentage of shrub coverage 0.0094 0.0041 0.0013 0.0175 

Percentage of herb coverage —0.0064 0.0028 —0.0118 —0.001 
Table 8 Comparison of morphological characteristics of species 

Spesies Body weight Claw length Culmen Wing length Body length Taillength Tarsus-metatarsus Finger length 
p (g) (mm) length (mm) (mm) (mm) (mm) length (mm) (mm) 
White-browed 
10.150 5.598 + 57.837 + 111.038 + 55.577 + 11.848 + 

Fulvetta (A.vini- 26 3 5.571+1.299° z b b b 24.946 +0.827° b 
pectus) 0.858 0.287 3.049 5.495 3.177 0.730 
eer 9.756+ 56294 54.8174 107.8734 52.397 + 11.618+ 
Fulvetta 63 6 680° 6.551 +1.386” 0 3210 24895 6 8310 3 444° 24.495 +0.931° 1.029? 
(A.cinereiceps) ： ： ` ： ： i 
Rusty-capped 
Fulvetta 31 15.803 + 6.861--0.884? 6.842 59.452 132.032 61.677 + 25.0793-0.961* 13.848 + 
(A.dubia 1.0867 0? 0.354? 1.748” 4.199° 3.250° 다 올가 1.403? 
ne Teese) 14.496-- 6751+  64208- 133000- 62.1924 11.927+ 
Fulvetta 120 1187: 5.828 +0.733° 6 362* men 3 043° 397a. 23.634 +0.830” Gave. 
(A.morrisonia) ： ` i ： ` i 
Golden-breasted 

Pa RS 7.126 a 47714 52.742+ 99645+ 48.9034 . — 10.1324 
Fulvetta 31 0.390* 4.716: 0.613 0.2364 1.751% 4.079 3.9443 22.792 +0.944 0.6614 
(A.chrysotis) : í i ： : i 
Rufous-winged 

RA 9.610+ a 5875+ 56.1004 104.100-  42.600+ 4 12.3004 
Fulvetta 10 0.985* 6.575 +0.834 0.196? 2.130 6.999% 2 836° 20.210+0.778 0.707° 
(A.castaneceps) ? ; ; ; : : 
Spectacled Fulv- 8.479 p; 9.6500 53.517 t 105.483 + 50.7414 명 10.635 + 
etta (A.ruficapilla) papal. 52911+1235 0.237* 2.0117 3.267* 3.081? ee eNOS 0.630° 
Same superscripts indicate non-significant differences. 

Approximately, 90.3% of the samples were correctly 


Linear discriminant analysis 

We conducted LDA to the eight morphological 
characteristics in Table 8, and got six linear discriminant 
functions. The cumulative proportion of LD1 and LD2 
reached 0.9284 (Table 9), indicating that they could account 
for 92.84% variability in the linear discriminant model. The 
proportions of the last four linear discriminant functions were 
small, indicating the weak influences they exert on the 
results. The absolute values of coefficients of body length 
(BL) and tarsus-metatarsus length (TML) in LD1 were high, 
indicating that LD1 primarily reflected the discriminant 
effects of body length and tarsus-metatarsus length. The 
absolute value of coefficient of wing length (WL) in LD2 was 
the largest coefficient, indicating that LD2 primarily reflected 
the discriminant effect of wing length (Table 10). 


Table9 Proportion of trace of linear discriminant analysis 


Functions LD1 LD2 LD3 


classified into different species via LDA. Samples were 
divided into four groupings and two of them were consisted of 
a single species, Golden-breasted Fulvettas and Rufous- 
winged Fulvettas, respectively. One group was consisted of 
Grey-cheeked Fulvettas and Rusty-capped Fulvettas, and the 
fourth group was consisted of White-browed Fulvettas, 
Streak-throated Fulvettas and spectacled Fulvettas (Table 11). 
The accuracy rates of LDA for each species were relatively 
high (91.20%-100%) except for White-browed Fulvettas 
(65.40%) and  Streak-throated — Fulvettas (74.60%). 
Misclassifications were found among White-browed Fulvettas, 
Streak-throated Fulvettas and spectacled Fulvettas, as well as 
between Grey-cheeked Fulvettas and Rusty-capped Fulvettas. 
Misclassification rate was relatively high between White- 
browed Fulvettas and Streak-throated Fulvettas (Table 11). 


LD4 LD5 LD6 





Proportion of trace 0.8523 0.0761 


0.0516 


0.0139 0.0051 0.001 
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Table 10 Coefficients of linear discriminant analysis 





Morphological characteristics LD1 LD2 LD3 LD4 LD5 LD6 

BW —0.8811885 0.3735774 —0.4774164 0.624646 0.1744586 0.2523426 
CUL —1.340453 0.1253922 1.2647371 —1.54594 1.4543998 —1.077508 
CLL 0.7308699 4.5448938 3.6643598 4.682241 6.6273355 5.5899434 
WL —14.8874514 —48.0531867 —12.1749091 30.530741 42.8527554 —2.7700184 
BL —22.2927648 —2.5958555 —15.5078565 4.240833 —26.0764156 15.6511549 
TL —3.5438221 9.4093755 —14.3753771 —29.584098 -7.0460453 —16.0740137 
TML 20.7689209 37.1068074 —933.2364027 16.419393 25.9016113 7.8645413 
FL 3.1495686 8.2214066 17.5987326 12.287153 —6.1706657 —17.0493372 


Table 11 Classification rates of linear discriminant analysis of morphological characteristics 


White-browed 
Fulvetta 
(A.vinipectus) 


Streak-throated 
Fulvetta 
(A.cinereiceps) 


spectacled 


Species Fulvetta 


White-browed Fulvetta 


65.409 
(A.vinipectus) à 


30.8096 3.8096 
Streak-throated Fulvetta 


i > 15.90% 74.60% 
(A.cinereiceps) 


9.50% 


Spectacled Fulvetta 
(A.ruficapilla) 


8.80% 91.20% 


Rusty-capped Fulvetta 
(A.dubia) 


Grey-cheeked Fulvetta 
(A.morrisonia) 


Golden-breasted Fulvetta 
(A.chrysotis) 


Rufous-winged  Fulvetta 
(A.castaneceps) 


DISCUSSION 


Spatial niche segregation in different scales 
According to the competitive exclusion principle, niche 
segregation was required among sympatric congeners to 
avoid competitive exclusion (Hardin, 1960; Levin, 1970; 
Schoener, 1974). Niche theory could well explain species 
coexistence in temperate forest (Nakashizuka, 2001). 
According to this study, the niche theory could also be used in 
subtropical montane forest to explain the coexistence of the 
seven fulvettas in Ailao Mountains, Ejia Twon, Yunnan 
Province. Through spatial niche segregation in different scales, 
fulvettas reduced inter-specific competition and promoted 
species coexistence. 

Habitat partitioning was relatively more important than other 
dimensions (Schoener, 1974). Habitat heterogeneity was vital 
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(A.ruficapilla) 


Rust d Grey-cheekea Coen Rufous-winged 
usty-cappe rey-cheeked sted ufous-winge 
Fulvetta Fulvetta Fulvetta 
(A.dubia) A.morrisonia Fund (A.castaneceps 
à M ) (A.chrysotis) i ps) 
93.50% 6.50% 
100% 


100% 


100% 


for habitat segregation (Vidus-Rosin et al, 2012). In this study, 
results of model selection based on GLM with Poisson error 
indicated that the elevation, shrub coverage and herb 
coverage were the key factors affecting fulvettas’ habitats 
selection (Table 6). Elevation analysis showed that niche 
segregation was detected through both mist-nets and point 
counts. No significant differences were found either among 
the three fulvettas occupying relatively high elevations or 
between Spectacled Fulvettas and Rusty-capped Fulvettas 
through mist-nets, whereas, through point counts, only one 
pair was found with non-significant differences. Chiang et al 
(2012) reported that elevation gradients might be the main 
factor in explaining the coexistence of species in spatial 
dimension. In the study area of this study, elevations were 
ranged from 800 m to 2 800 m and both the climate and forest 
types were showed with obvious vertical variations, therefore, 
the habitat heterogeneity along elevation gradients offered 


opportunities of these fulvettas to choose different habitats. 

Niche overlaps were also existence for some fulvettas 
even though segregation in elevation had been detected. 
Schoener (1983) found that when species were similar in 
one dimension, resource differentiation would occur in other 
dimensions to reduce inter-specific competition. Vertical 
height was an important component of spatial niche. Studies 
on tits had shown that segregation in vertical height could 
facilitate species' coexistence (Song, 1983; Yang et al, 
2012). Through analysis of roost site height, we found niche 
overlaps among  White-browed  Fulvettas, Spectacled 
Fulvettas and  Streak-throated — Fulvettas. However, 
Spectacled Fulvettas had significant differences with the 
other two in elevation. Golden-breasted Fulvettas was 
widely overlapped with White-browed Fulvettas and Streak- 
throated Fulvettas in elevation documented by mist-nets. 
They had significant differences in roost site height with the 
other two species. Similar results were also found between 
Rufous-winged Fulvettas and Streak-throated Fulvettas 
documented by point counts. Vertical height reflects the 
activity space chosen by birds. Because of their unique 
physiological and activity pattern, birds have high demand 
for energy supply. Foraging behavior accordingly occupies a 
large proportion in birds’ activities. Therefore, differentiation 
in roost site height mainly reflects segregation in foraging 
height. Studies on tits (Liu et al, 1989), prinias (Zhou & Fang, 
2000), hummingbirds (Lara et al, 2011) and other bird 
communities (Gao & Yang, 1991) found that foraging height 
segregation reduced inter-specific competition and 
facilitated coexistence. Hence, we assumed that foraging 
height segregation permitted fulvettas that widely 
overlapped in elevation to relax the intensity of inter-specific 
competitive interactions. 

Through analysis of elevation and vertical height, we did 
not found significant differences between White-browed 
Fulvettas and Streak-throated Fulvettas documented by 
mist-nets. However, significant differences were then 
detected in the key factors of shrub coverage and herb 
coverage (Table 5). Vegetation coverage was an important 
component of  microhabitat chosen by animals. 
Segregation in  microhabitat could facilitate species’ 
coexistence (Dammhahn et al, 2013; Traba et al, 2013). 
Study on Blue Eared Pheasant (Crossoptilon auritum) 
found that concealment condition provided by vegetation 
coverage had a significant influence on birds' habitat 
selection (Liu et al, 2005). Study of breeding ecology 
showed that fulvettas' nests were primarily located in 
undergrowth consisted of shrubs and herbs (Lee et al, 
2010; Gong, 1994; Huang et al, 1988; Zhou, 1989). 
Therefore, segregation in this scale allowed fulvettas to 
use different concealment conditions in their home range 
and avoided complete niche overlap accordingly. 


Morphological differentiations 

Morphological characteristics provide an insight into the 
ecology of animals (Landmann & Winding, 1993; Miles et 
al, 1987). It reflects the adaptions to the environmental 


conditions consisted of abiotic and biotic factors during 
animals' life history (Martin, 2001). We detected significant 
differences in eight morphological characteristics of the 
seven fulvettas (Table 8) and four qualitative groupings 
were resulted from LDA (Table 11). Numerous of studies 
on the relationships between birds’ niche use and 
morphology found that differences in morphological 
characteristics led to differentiation in their competitive 
abilities (Kalinoski, 1975; Gao et al, 1997), their foraging 
behaviors (Salewski et al, 2003), their resource 
preferences (Hill & Lein, 1988) and the varieties, sites and 
vertical height of trees they chose (Alatalo, 1981; Richards 
et al, 2000; Salewski et al, 2003). In this study, LDA 
primarily reflected discriminant effects of body length, 
tarsus-metatarsus length and wing length. These three 
characteristics have great influences on bird's activities. 
Studies on herons (Wen et al, 1998; Ye et al, 2006; Zhu et 
al, 1998), hummingbirds (Lara et al, 2011) and 
woodpeckers (Gao et al, 1997) showed that body size 
played a key role in deciding birds' competitive ability. 
Wing length decides birds' flying ability. Tarsus-metatarsus 
length has a strong correlation with birds' behavior and 
influences habitat selection (Liu et al, 2013). Therefore, we 
predicted that fulvettas of different groupings present 
different features in resource use which helped them to 
realize niche partition. Misclassifications within groups 
indicated the similarity in resource use among group 
members. However, we also detected segregations among 
fulvettas within each group. For example, in the first group, 
White-browed Fulvettas was different from Streak-throated 
Fulvettas in microhabitat use, and they had significant 
differences with Spectacled Fulvettas in elevation. 
Significant differences were also detected in elevation and 
roost site height between Grey-cheeked Fulvettas and 
Rusty-capped Fulvettas of the second group (Table 11). 

In conclusion, different morphological characteristics of 
each group benefited fulvettas to partition resource in the 
overlapped regions and spatial niche segregation relaxed the 
intensity of inter-specific interactions among members within 
each group. However, the misclassifications indicated that 
intense competition could still occur among members within 
each groups in the overlapped regions and these 
competitions might become new selective pressures 
facilitating further differentiations. 

The findings of this study showed that niche theory could 
explain the coexistence mechanisms of the seven fulvettas 
in the subtropical montane, Ejia Town, Yunnan Province. 
Through elevation, roost site height and vegetation coverage 
partitioning, these seven sympatric fulvettas realized spatial 
niche segregation in different scales. Combined with 
differentiations in resource use due to different 
morphological characteristics, they were able to minimize 
the intensity of inter-specific interactions and promote the 
coexistence. Moreover, the mutability and unpredictability of 
environment, the interference effects, the migration of 
species and the heterogeneity provided by large 
environmental gradient might also play important roles and 
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should be considered along with the competition effects 
when exploring the coexistence mechanisms of sympatric 
congeners. 
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ABSTRACT 


In reptiles, habitat selection is the process whereby 
suitable habitat is selected that optimizes 
physiological functions and behavioral performance. 
Here, we used the brown forest skink (Sphenomor- 
phus indicus) as a model animal and examined 
whether the frequency of active individuals, envi- 
ronmental temperature, illumination of activity area, 
and habitat type vary with different age classes. We 
surveyed the number of active individuals and 
measured environmental variables at Baiyunshan 
Mountain in Lishui, Zhejiang, China. We found no 
difference in the activity frequency of adult and juve- 
nile S. indicus; the activity pattern of active individu- 
als was bimodal. The mean environmental 
temperature selected by adults was higher than that 
selected by juveniles. The environmental 
temperature of active areas measured at 0900-1000 
h and 1100-1200h was higher than at 1400-1500h; 
illumination of the active area at 1000-1200h was 
also higher than at 1400h-1600 h. The number of 
active individuals, the environmental temperature 
and illumination of activity areas showed 
pairwise positive correlation. There was a difference 
in habitat type between juveniles and adults whereby 
juveniles prefer rock habitats. We predict that active 
S. indicus select optimal habitats with different 
environmental temperatures and types to reach the 
physiological needs particular to their age classes. 


Keywords: Sphenomorphus indicus; Age-related 
selection difference; Habitat selection; Temperature; 
Illumination; Habitat type; Activity frequency 


INTRODUCTION 


Habitat selection is the process of animals choosing optimal 
habitat (Buckland et al, 2014; Johnson, 1980). Studies on 
habitat selection can help us understand animal behavior, 
population dynamics, propagation, survival and wildlife 
conservation (Buckland et al, 2014; Manly et al, 2002; 
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Strickland & McDonald, 2006). Because habitat selection 
variation occurs with spatial dimensional change, we can 
acquire information about habitat selection by examining 
animal distributions at specific spatial scales (De La Cruz et al, 
2014; Hódar et al, 2000; Oppel et al, 2004). Studies on habitat 
selection have mainly focused on the three levels of popula- 
tion, home range and microhabitat (Thomas & Taylor, 1990). 
However, studies at different levels are important in 
determining the factors vital to habitat utilization and effective 
conversation measures (Beasley et al, 2007; Buckland et al, 
2014; Razgour et al, 2011). 

Habitat selection in most ectotherms is determined by spe- 
cific thermal preferences, social interactions, predator avoid- 
ance and other factors (Angilletta et al, 2002; Dochtermann et 
al, 2012; Downes & Shine, 1998; Hall et al, 1992). Ecto- 
therms are able to optimize their physiological and functional 
performance through habitat selection, and studies at the 
microhabitat level are important in understanding the mecha- 
nisms of selection. As a representative species and ectother- 
mic group, reptiles not only need to avoid predators during 
activity, but more importantly, their behavior and physiological 
capabilities are significantly affected by temperature. It is vital 
that reptiles find suitable heterogeneous thermal environments 
to improve fitness (Angilletta et al, 2002; Castilla et al, 1999; 
Huey, 1982; Webb, 1996). During microhabitat selction, the 
activities of a reptile are under the influence of other repetiles or 
other species within the environment (Manicom & Schwarzkopf, 
2011; Stamps & Tanaka, 1981; Wasko & Sasa, 2012). 
Sometimes, competition over environmental resources occurs 
(e.g., food, thermal resources) (Amarasekare & Coutinho, 2014; 
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Marler et al, 1995; Hart & Marshall, 2012) and differences in 
intraspecies or interspecies habitat selection will follow. 

The brown forest skink (Sphenomorphus indicus) is a member 
of the family Scincidae and mainly distributed across sourthen 
China in dark and wet grass, rock piles or cracks in cliffs (Huang, 
1999). They are active away from their caves from early April to 
late October, and are more likely to be observed in shady areas 
during the morning and afternoon, and occasionally at noon. 
Here, we used line transects to investigate the quantity and type 
of habitat used by S. indicus in hilly road areas of Baiyunshan 
Mountain, Lishui, China. Environmental temperatures and illumi- 
nation intensity were measured to explore characteristics and 
underlying mechanisms of habitat selection by S. indicus of 
different ages. Our results will aid the conservation of wild 
populations of this species. 


MATERIALS AND METHODS 


Experimental design 

From May to early July 2014 we used line transect methods to 
investigate S. indicus and their habitat at Baiyunshan Mountain, 
Lishui, Zhejiang, China (N27°49'17.50", E120*32'14.89"). The 
full length of the line transect was 4 km, and the altitude at the 
start and end points was 140 m and 360 m, respectively. During 
the survey observers proceeded at an even pace along hilly 
road. Active individuals sighted within 2 m of the road along 
the cliff (width=2 m) and environmental parameters of habitat 
were recorded. Seventeen surveys were randomly carried out 
during 0900h-1230h (10 surveys) and 1230h-1600h (7 
surveys) over 6 days. The weather and temperature during 
surveys are shown in Table 1. Active individuals were 
determined as adults (SVL> 65 mm) or juveniles (SVL«65 mm) 
(Ji & Du, 2000) by visually estimating body size. The snout-vent 
length (SVL) of randomly captured individuals was measured 


Table1 Weather conditions during field surveys 


Date Time Weather Temperature (9C) 





16/05 1230h-1600h Sunny 22.8 (20.7-26.5) 


21/05 0900h-1230h Sunny 24.1 (17.1-31.8) 
24/05 0900h-1230h Cloudy 24.8 (20.4-31) 
25/05 0900h-1230h Sunny 26.4 (23.2-32) 

27/05 0900h-1600h Overcast to cloudy 24.1 (22.3-29.4) 
01/06 0900h-1600h Sunny 26.6 (22.7-32.3) 
07/06 0900h-1600h Overcast to cloudy 26.9 (22.0-34.3) 
25/06 0900h-1600h Overcast to sunny 25.4 (20.8-32) 

06/07 0900h-1230h Showers to cloudy 29.6 (26.0-34.0) 
07/07 0900h-1600h Showers to cloudy 26.4 (26.0-36.0) 


09/07 0900h-1600h Cloudy to sunny 29.6 (26.0-36.0) 
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by digital square caliper (0.01 mm, Mitutoyo, Japan) and the 
accuracy of visual estimation was verified. Before being re- 
leased to the original capture area, the 1*' toe on the left hind 
legs was clipped in individuals to avoid repeated sampling. 
Because sex can not be determined by visual estimation, 
gender-related differences were not considered here. The 
habitat of observed S. indicus was determined whereby grass 
or leaf litter layer was categorized as grass habitat and bare 
rock or ditch was categorized as rock habitat. An infrared 
radiation thermometer (UT301A, Uni-T, China) was used to 
measure the temperature of the spot where a skink was 
located, defined by a circle with a 10 cm diameter around the 
animal; the environmental temperature of this habitat was the 
average value of 10 measurements. A sensitive probe was 
planted at the center of this circle to measure illumination 
intensity (lux). 


Statistical analysis 

Data were analyzed using STATISTICA (v6.0). The normality 
and homogeneity of variance of all data were checked by 
Kolmogorov-Simirnov tests and Bartlet tests, respectively; no 
parameters required transformation before statistical analyses. 
The percentage of active S. indicus during a survey was equal 
to the numbers of observed active S. indicus during a survey 
divided by the total numbers of observed S. indicus during a 
survey, multiplied by 100. Contingency tables and G-tests 
were used to analyze age-related activity patterns and habitat 
type, and activity frequencies during each time period. Two- 
factor ANOVA was used to determine the effects of age and 
survey time on environmental temperature and illumination 
intensity of the active area. One-factor ANOVA was used to 
determine the effects of age-related individual differences and 
habitat type on environmental temperature and illumination 
intensity of the active area. Partial correlation was used to 
check correlations between numbers of active individuals 
during a survey period, environmental temperature and 
illumination intensity of the active area. Group differences were 
compared using Tukey's post hoc test. All data were expressed 
as mean+SE and statistical significance was taken at a=0.05. 


RESULTS 


We collected 404 S. indicus (211 juveniles and 193 adults), 
and acquired parameters for habitat type, environmental 
temperature and illumination intensity in active areas. Body 
size was measured in 96 individuals (31 juveniles and 65 
adults). The SVL of adults (74.7+ 0.6 mm, 65.8-91.3 mm) was 
longer than that of juveniles (60.5+0.7 mm, 46.0-64.7 mm) (F4, 
947205.82, P«0.001) and no overlap was found between the 
SVL ranges of adults and juveniles. The activity patterns of 
adults and juveniles during each survey period were 
comparable (G=9.04, df-6, P=0.171) (Figure 1A). However, 
when data from adults and juveniles were grouped together, 
differences in the proportions of active individuals during each 
survey period were found (G=178.22, df-6, P«0.001). For 
example, two activity peaks occurred at 1000-1200 h and 
1300-1500 h, and individuals were more active during 1000- 
1200 h (Figure 1B). 
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Figure 1 Activity frequency of juvenile and adult Sphenomorphus 
indicus (A) and proportion of active individuals at survey times (B) 
Numbers above bars indicate sample size for each time. 


Both age and survey time influenced the environmental 
temperature of areas of activity in S. indicus, but the age x 
urvey time interaction did not. The environmental temperature 
of areas of adult activity (24.8+0.2 °C) was higher than for 
juveniles (24.140.2 °C) (Tukey's test, P<0.001). The 
environmental temperatures during 0900-1000 h and 1100- 
1200 h were both higher that during 1400-1500 h (Tukey’s 
test, both P<0.03); no differences were found across other 
survey periods (Tukey's test, all P>0.05). 

Survey time influenced illumination intensity of the areas 
frequented by active S. indicus, but age and the agex survey 
time interaction did not. The illumination intensity during 1000- 
1200 h was higher than during 1400-1600 h (Tukey’s test, all 
P<0.04); no differences were found among other survey 
periods (Tukey’s test, all P>0.05) (Figure 2, Table 2). During 
each survey period, the number of active individuals, envi- 
ronmental temperature and illumination intensity of the active 
area were positively correlated (all P«0.04). 


Table2 Statistical results of environmental temperature and 
illumination in areas where Sphenomorphus indicus are active 








Environmental tempera- AE : 
Illumination intensity 
ture 
E df P F df P 
Age 3.88 1,390 <0.001 1.10 1,390 0.295 
Measurement time 8.00 6, 390 <0.01 4.58 6,390 <0.001 
Interaction 0.92 6,390 0.480 0.97 6,390 0.448 
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Figure 2 Environmental temperature (A) and illumination inten- 


sity (B) of areas where Sphenomorphus indicus are active over 
different survey times 


Habitat type in the activity areas did not affect either 
environmental temperature (F1402=0.05, P=0.829) or 
illumination intensity (F1402=2.49, P=0.115). Habitat selection 
between adults and juveniles was different (G=4.70, df=1, 
P<0.04) (Figure 3). 
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Figure 3 Proportion of habitat types selected by active juvenile 
and adult Sphenomorphus indicus 


DISCUSSION 


The activity patterns of reptiles is affected by many biotic and 
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abiotic factors (Cloudsley-Thompson, 1961; Oishi et al, 2004), 
including age-related energy requirements among different 
individuals (Liu et al, 2008; Mautz & Nagy, 1987; Oishi et al, 
2004; Paulissen, 1987; ). For example, the daily activity of 
common lizard (Lacerta vivipara) in juveniles is started later 
and shorter than in adults (Liu et al, 2008). Here, activity 
patterns were similar between adult and juvenile S. indicus 
(Figure 1A), indicating that because the habitat of S. indicus 
is relatively shady, wet and closed with high food abundance, 
skinks are less likely to compete with each other over food. It 
appears that niche separation over activity time distribution is 
not induced among individuals of different ages in this species 
and this is consistent with other studies on diurnal activity 
rhythms in other diurnal lizards (Liang et al, 2006; Liu et al, 
2008; Wen & Zou, 2002;). Activity rhythms with double peaks 
were found in S. indicus during summer. The highest 
percentage of active individuals occurred during 1000-1200 h, 
and this decreased during 1200-1300 h and increased again 
during 1300-1500 h (Figure 1B). Differences between our 
results and those reported by Wang (1964) on S. indicus in 
Hangzhou are probably the result of local weather and geo- 
graphic conditions. 

As typical ectotherms, the biological function and behavioral 
performance of reptiles are significantly influenced by the 
thermal environment (Angilletta et al, 2002; Angilletta, 2009; 
Shu et al, 2010; Vermunt et al, 2014). Extreme temperatures 
can be harmful to individuals and good fitness requires a 
suitable temperature range (Clusella-Trullas & Chown, 2014; 
Shen et al, 2013). We found that the mean environmental 
temperature of the active areas frequented by adult S. indicus 
was 0.7 °C higher than the areas used by juveniles (Figure 2); 
a similar phenomenon has been observed in other thermal 
biology studies on the Iberian wall lizard (Podarcis hispanica 
atrata), the Mongolia racerunner (Eremias argus) and the 
multi-ocellated racerunner (Eremias multiocellata) (Castilla & 
Bauwens, 1991; Tang et al, 2013; Xu & Ji, 2006). These age- 
related differences in body temperature and environmental 
temperature of habitat may be correlated with body size be- 
cause a small body size is usually characterized by higher 
thermo-conductivity but weaker heat storage capacity (Castilla 
& Bauwens, 1991; Tang et al, 2013; Xu & Ji, 2006). Moreover, 
the differences in feeding activity and predator avoidance 
during development may induce variation in thermal environ- 
ments (Christian & Bedford, 1995; Hertz et al., 1993). The 
similar illumination intensities among habitats of individuals at 
different ages (Figure 2) indicate that illumination intensity 
does not affect habitat selection. However, temporal rhythms 
were found in both adult and juvenile S. indicus, i.e., the envi- 
ronmental temperature and illumination intensity of the active 
areas were both higher in mornings (0900-1200 h) than after- 
noons (1400-1600 h) (Figure 2), indicating that the rhythms of 
S. indicus choosing environmental factors are consistent with 
diurnal activity patterns. Positive correlations among the num- 
ber of active individuals, environmental temperature and 
illumination intensity also indicate that the activity rhythms of 
S. indicus are interwoven with environmental thermal and 
illuminative factors. 


32 WWW.ZOOres.ac.cn 


Reptiles choose optimal habitats to obtain more thermal re- 
sources and a better chance of obtaining food (Wang, 1964; 
Strickland & McDonald, 2006; Buckland et al, 2014). Our 
study found that juvenile S. indicus preferred rock habitats 
(6496) during diurnal activity (Figure 3), consistent with Wang 
(1964). However, neither the environmental temperature nor 
illumination intensity differed between rock habitat and grass 
habitat, suggesting that the reason juveniles choose rock 
habitat is not based on acquiring more heat and illumination. 
The optimal physiological function and behavioral 
performance of lizards changes during development (Tang et 
al, 2013; Xu & Ji, 2006), therefore, lizards at different ages 
should choose different habitats to meet their physiological 
requirements and avoid predators (Irschick et al, 2000). For 
example, two species of anole lizards (Anolis lineatopus and 
Anolis gundlachi) in India choose habitats located in low and 
narrow areas to avoid predators (Irschick et al, 2000). Be- 
cause the dark brown body color of S. indicus (Huang, 1999) 
is close to the color of rocks, we assume that juveniles prefer 
moving on rocks so as to avoid predators. We therefore 
predict that active S. indicus select optimal habitats with 
different environmental temperatures and types to reach the 
physiological needs particular to their age classes. 


ACKNOWLEDGEMENTS 


We are grateful to Shi-Liang QIN, Yu YU, Xiao-Wei WANG, Jian-Yang XIE, 
Li-Shuang ZHANG, Qiong-Lu ZHANG, Xue-Ya ZHANG, Xiao-Ning ZHOU 
for their help during this research. 


REFERENCES 


Amarasekare P, Coutinho RM. 2014. Effects of temperature on intraspecific 
competition in ectotherms. American Naturalist, 184(3): E50-E65. 

Angilletta MJ Jr. 2009. Thermal Adaptation: A Theoretical and Empirical 
Synthesis. Oxford: Oxford University Press. 

Angilletta MJ Jr, Niewiarowski PH, Navas CA. 2002. The evolution of 
thermal physiology in ectotherms. Journal of Thermal Biology, 27(4): 249- 
268. 

Beasley JC, Devault TL, Retamosa MI, Rhodes OE. 2007. A hierarchical 
analysis of habitat selection by raccoons in northern Indiana. Journal of 
Wildlife Management, 71(4): 1125-1133. 

Buckland S, Cole N C, Godsall B, Rodriguez-Pérez JR, Gallagher LE, 
Henshaw SM, Harris S. 2014. Habitat selection of the Mauritian lowland 
forest day gecko at multiple spatial scales: a baseline for translocation. 
Global Ecology and 
10.1016/j.gecco.2014.06.001. 


Conservation, 1: 71-79, doi: 


Castilla AM, Bauwens D. 1991. Thermal Biology, microhabitat selection, 
and conservation of the insular lizard Podarcis hispanica atrata. Oecologia, 
85(3): 366-374. 

Castilla AM, Van Damme R, Bauwens D. 1999. Field body temperatures, 
mechanisms of thermoregulation and evolution of thermal characteristics in 
Lacertid lizards. Natura Croatica, 8(3): 253-274. 

Christian KA, Bedford GS. 1995. Seasonal changes in thermoregulation by 


the frillneck lizard, Chlamydosaurus kingii, in tropical Australia. Ecology, 
76(1): 124-132. 


Cloudsley-Thompson JL. 1961. Rhythmic Activity in Animal Physiology and 
Behaviour. New York: Academic Press. 

Clusella-Trullas S, Chown SL. 2014. Lizard thermal trait variation at multi- 
ple scales: a review. Journal of Comparative Physiology B, 184(1): 5-21. 


De La Cruz SEW, Eadie JM, Miles AK, Yee J, Spragens KA, Palm EC, 
Takekawa JY. 2014. Resource selection and space use by sea ducks 
during the non-breeding season: implications for habitat conservation 
planning in urbanized estuaries. Biological Conservation, 169: 68-78. 


Dochtermann NA, Jenkins SH, Swartz MJ, Hargett AC. 2012. The roles of 
competition and environmental heterogeneity in the maintenance of behav- 
ioral variation and covariation. Ecology, 93(6): 1330-1339. 


Downes 5, Shine R. 1998. Heat, safety or solitude? Using habitat selection 
experiments to identify a lizard’s priorities. Animal Behaviour, 55(5): 1387- 
1396. 


Hall CAS, Stanford JA, Hauer FR. 1992. The distribution and abundance of 
organisms as a consequence of energy balances along multiple environ- 
mental gradients. Oikos, 65(3): 377-390. 


Hart SP, Marshall DJ. 2012. Advantages and disadvantages of interference- 
competitive ability and resource-use efficiency when invading established 
communities. Oikos, 121(3): 396-402. 
Hertz PE, Huey RB, Stevenson R. 1993. Evaluating temperature regulation 
by field-active ectotherms: the fallacy of the inappropriate question. The 
American Naturalist, 142(5): 796-818. 


Hódar JA, Pleguezuelos JM, Poveda JC. 2000. Habitat selection of the 
common chameleon (Chamaeleo chamaeleon) (L.) in an area under devel- 
opment in southern Spain: implications for conservation. Biological Conser- 
vation, 94(1): 63-68. 

Huang QY. 1999. Sphenomorphus indicus (Gray, 1853). In: Zhao EM, Zhao 
KT, Zhou KY. Fauna Sinica, Reptile Vol. 2, Squamata, Lacertilia. Beijing: 
Science Press, 340-349. (in Chinese) 

Huey RB. 1982. Temperature, physiology, and the ecology of reptiles. /n: 
Gans C C, Pough FH. Biology of the Reptilia, Vol. 12. New York: Academic 
Press, 25-91. 

Irschick DJ, Macrini TE, Koruba S, Forman J. 2000. Ontogenetic differ- 
ences in morphology, habitat use, behavior, and sprinting capacity in two 
West Indian Anolis lizards. Journal of Herpetology, 34(3): 444-451. 

Ji X, Du WG. 2000. Sexual dimorphism in body size and head size and 
female reproduction in a viviparous skink, Sphenomorphus indicus. Zoo- 
logical Research, 21(5): 349-354. (in Chinese) 

Johnson DH. 1980. The comparison of usage and availability measur- 
ements for evaluating resource preference. Ecology, 61(1): 65-71. 

Liang WB, Zhang YX, Su P, Long Q, Huang JQ. 2006. Observation on time 
budget of Shinisaurus crocodilurus in captivity. Sichuan Journal of Zoology, 
25(2): 264-266. (in Chinese) 

Liu P, Liu ZT, Li DW, Zhao WG. 2008. Diurnal activity rhythm and time 
budget of Lacerta vivipara in simulated habitat. Chinese Journal of Ecology, 
27(12): 2146-2152. (in Chinese) 

Manicom C, Schwarzkopf L. 2011. Diet and prey selection of sympatric 
tropical skinks. Austral Ecology, 36(5): 485-496. 

Manly BFJ, McDonald LL, Thomas DL, McDonald TL, Erickson WP. 2002. 
Resource Selection by Animals: Statistical Design and Analysis for Field 


Studies. 2"" ed. Netherlands: Kluwer Academic Publishers. 


Marler CA, Walsberg G, White ML, Moore M. 1995. Increased energy 


expenditure due to increased territorial defense in male lizards after pheno- 
typic manipulation. Behavioral Ecology and Sociobiology, 37 (4): 225-231. 


Mautz WJ, Nagy KA. 1987. Ontogenetic changes in diet, field metabolic 
rate, and water flux in the herbivorous lizard Dipsosaurus dorsalis. Physio- 
logical Zoology, 60(6): 640-658. 

Oishi T, Nagai K, Harada Y, Naruse M, Ohtani M, Kawano E, Tamotsu S. 
2004. Circadian rhythms in amphibians and reptiles: ecological implications. 
Biological Rhythm Research, 35(1-2): 105-120. 


Oppel S, Schaefer HM, Schmidt V, Schróder B. 2004. Habitat selection by 
the pale-headed brush-finch (Atlapetes pallidiceps) in southern Ecuador: 
implications for conservation. Biological Conservation, 118(1): 33-40. 


Paulissen MA. 1987. Optimal foraging and intraspecific diet differences in 
the lizard Cnemidophorus sexlineatus. Oecologia, 71(3): 439-446. 


Razgour O, Hanmer J, Jones G. 2011. Using multi-scale modelling to 
predict habitat suitability for species of conservation concern: the grey long- 
eared bat as a case study. Biological Conservation, 144(12): 2922-2930. 


Shen JW, Meng FW, Zhang YP, Du WG. 2013. Field body temperature and 
thermal preference of the big-headed turtle Platysternon megacephalum. 
Current Zoology, 59(5): 626-632. 

Shu L, Zhang QL, Qu YF, Ji X. 2010. Thermal tolerance, selected body 
temperature and thermal dependence of food assimilation and locomotor 
performance in the Qinghai toad headed lizard, Phrynocephalus viangalii. 
Acta Ecologica Sinica, 30(8): 2036-2042. (in Chinese) 

Stamps JA, Tanaka S. 1981. The relationship between food and social 
behavior in juvenile lizards. Copeia, 1981(2): 422-434. 

Strickland MD, McDonald LL. 2006. Introduction to the special section on 
resource selection. Journal of Wildlife Management, 70(2): 321-323. 


Tang XL, Yue F, He JZ, Wang NB, Ma M, Mo JR, Chen Q. 2013. Ontoge- 
netic and sexual differences of thermal biology and locomotor performance 
in a lacertid lizard, Eremias multiocellata. Zoology, 116(6): 331-335. 


Thomas DL, Taylor EJ. 1990. Study designs and tests for comparing resource 
use and availability. Journal of Wildlife Management, 54(2): 322-330. 


Vermunt A, Hare KM, Besson AA. 2014. Unusual change in activity pattern 
at cool temperature in a reptile (Sphenodon punctatus). Journal of Thermal 
Biology, 42: 40-45. 


Wang PC. 1964. Ecology of four lizards in Hangzhou (1): distribution, 
activity rhythm and food habits. Chinese Journal of Zoology, 6(2): 70-76. (in 
Chinese) 


Wasko DK, Sasa M. 2012. Food resources influence spatial ecology, 
habitat selection, and foraging behavior in an ambush-hunting snake 
(Viperidae: Bothrops asper): an experimental study. Zoology, 115(3): 179- 
187. 


Webb J. 1996. Ecology of A Threatened Snake Species, Hoplocephalus 
bungaroides (Elapidae). Ph. D. thesis, University of Sydney, Sydney. 


Wen CY, Zou PZ. 2002. Preliminary studies on population ecology and 
diet of Eumeces chinensis in the north of Guangdong province. Journal 
of Guangzhou University(Natural Science Edition), 1(3): 19-22. (in 
Chinese) 


Xu XF, Ji X. 2006. Ontogenetic shifts in thermal tolerance, selected body 
temperature and thermal dependence of food assimilation and locomotor 
performance in a lacertid lizard, Eremias brenchleyi. Comparative Biochem- 
istry and Physiology Part A: Molecular & Integrative Physiology, 143(1): 
118-124. 


Zoological Research 36(1): 29-33, 2015 33 


ZOOLOGICAL RESEARCH 


Microsatellite analysis of genetic diversity and population 
structure of freshwater mussel (Lamprotula leai) 


Jin-Jin MIN', Rong-Hui YE?*, Gen-Fang ZHANG?, Rong-Quan ZHENG” 


1 Institute of Ecology, Zhejiang Normal University, Jinhua Zhejiang 321004, China 


? Jinhua Polytechnic, Jinhua Zhejiang 321004, China 


ABSTRACT 


Lamprotula leai is one of the most commercially 
important freshwater pearl mussels in China, but 
there is limited data on its genetic diversity and 
population structure. In the present study, 119 
individuals from four major geographical populations 
were investigated using 15 microsatellite loci 
identified via cross-species amplification. A total of 
114 alleles were detected, with an average of 7.6 
alleles per locus (range: 2 to 21). Among the four 
stocks, those from Hung-tse Lake and Poyang Lake 
had the lowest (0.412) and highest (0.455) observed 
heterozygosity respectively. The polymorphism 
information content (PIC) ranged from 0.374 to 
0.927 (mean: 0.907). AMOVA showed that 12.56% 
and 44.68% genetic variances were among 
populations and within individuals, respectively. 
Pairwise Fst ranged from 0.073 to 0.146, indicating 
medium genetic differentiation among the populations. 
In aggregate, our results suggest that inbreeding is a 
crucial factor accounting for deviations from Hardy— 
Weinberg equilibrium at 12 loci. Moreover, the 
genetic distance among four stocks ranged from 
0.192 to 0.890. Poyang Lake and Hung-tse Lake 
were clustered together, joined with Dongting Lake 
and Anqing Lake. Given that specimens from Hung- 
tse Lake showed the highest average allele richness, 
expected heterozygosity and PIC, this location may 
be the source of the highest quality germplasm 
resources and the stock from this area may be the 
best for future breeding efforts. 


Keywords: Lamprotula leai; Freshwater mussel; 
Genetic diversity; Population structure; Microsatellite loci 


INTRODUCTION 


Genetic diversity is important for sustainable exploitation of 
cultured resources (Afanas’ev et al, 2006; Laikre et al, 2005), 
especially as the exploitation of aquatic stocks and 
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environmental degradation of habitat becomes more 
commonplace. Higher levels of genetic diversity among these 
stocks often grants them greater ability to respond to 
environmental changes, artificial selection and pathogen 
infection, all of which tend to occur in during intensified 
aquaculture (Liu & Yao, 2013; Wu et al, 2013). For example, 
Lamprotula leai is an endemic species distributed in large and 
medium rivers and lakes across China (Hu, 2005) that are 
widely used in pearl aquaculture and indigenous handicrafts 
due to their large shell, strong ability to secrete pearls, and 
thick nacre (Wang et al, 2007). Likewise, this mussel is widely 
consumed as food throughout China (Liu et al, 1979). Despite 
the importance of this species, little is known about the current 
state of its genetic diversity, except that among cultured species 
there is germplasm degradation and general declines genetic 
diversity, potentially due to intensive farming and the method of 
cultivation (Ling, 2005). Due to the intensive cultivation of this 
species, it is now listed as a first-class protected aquatic wildlife 
species in Anhui province and second-class protected species 
in Hubei province (Xu et al, 2012). 

To date, most of the research into L. leai has focused on 
age and growth (Ling et al, 2005), conservation biology (Ling, 
2005), embryonic development (Zhang et al, 2009), abnormal 
development and effective accumulated temperature of 
parasitic glochidium (Zhang et al, 2010a), and growth and 
development of juvenile mussels (Zhang et al, 2010b). 
However, little work has been done on this species’ genetic 
diversity and population structure. One strategy to improve 
our understanding of population structure and genetic 
diversity of many aquaculture species is molecular analysis 
(Liu & Corde, 2004). Microsatellite markers have been shown 
to be suitable tools to assess genetic diversity because of 
their intrinsic genetic characteristics, including high 





Received: 28 May 2014; Accepted: 03 December 2014 

Foundation items: Public Welfare Projects of Technology Office in 
Zhejiang Province (2011C32SA700049); Major Science and Techno- 
logy Specific Projects of Zhejiang Province (2012C12907-9); Science 
and Technology Plan Projects of Jinhua City (2011A22020) 
“Corresponding authors, E-mails: yeronghui@163.com; zhengrq@zjnu.cn 


DOI:10.13918/j.issn.2095-8137.2015.1.34 


Zoological Research 36(1): 34-40, 2015 


polymorphism, stability and specificity, and co-dominant 
inheritance (Sun et al, 2008). More promisingly for studying L. 
leai, microsatellites have been previously used to analyze 
genetic diversity in several mollusks, including Hyriopsis 
cumingii (Ji, 2007), Cristaria plicata (Ji et al, 2007), Pinctada 
martensii (Yan et al, 2009), and Anodonta woodiana (Wang et 
al, 2011a). To date, 18 pairs of microsatellite primers have 
been isolated using magnetic bead hybridization and 5'- 
anchored PCR methods (Xu et al, 2011, 2012). In this study, 
we used these microsatellites to analyze the genetic 
relationships among four different stocks of L. leai and provide 
a novel theoretical basis for genetic resource protection and 
genetic management. 


MATERIALS AND METHODS 


Sample collection and DNA isolation 

Freshwater mussel specimens originating from four 
geographical locations across China—Poyang Lake (PY), 
Dongting Lake (DT), Hung-tse Lake (HZ), and Anqing Lake 
(AQ)—were obtained from the Weiwang Pearl Cultivation 
Base in Zhejiang Province (Figure 1), with 31, 31, 28, and 29 
samples respectively from PY, DT, HZ and AQ (Hale et al, 
2012). Genomic DNA was extracted from muscular tissue 
following a standard phenol: chloroform protocol as published 
previously (Sambrook & Russell, 2001). 













Dongting Lake 


Poyang Eake 














Figure 1 Sampling locations of four Lamprotula leai stocks in China 


Microsatellite amplification in L. leai 

We used microsatellite primers published for H. cumingii (Bai 
et al, 2009; Li et al, 2007; Luo, 2006; Wang et al, 2006; Xu et 
al, 2010; Zhu et al, 2010) and other species closely related to 
L. leai (Ji, 2007; Launey & Hedgecock, 2001). A total of 107 
candidate primer pairs (see supplemental Table 1, supporting 
information of http://www.zoores.ac.cn/) were synthesized 
(Shanghai Sangon Company) and each microsatellite was 
amplified in a 25 uL PCR containing 50 ng of DNA, 1 uL each 
of 10 umol/L primer, 2.5 uL of 10x buffer, 2 uL of dNTP 
(10 mmol/L), 1 U of Taq polymerase (5 U/uL), and 17.3 uL of 
ddH2O. PCR was conducted under the following conditions: 4 
min denaturation at 94 °C 32 cycles of 30 s at 94 °C, 30 s at 
specific annealing temperatures, and 30s at 72°C; and a 
final extension at 72°C for 10 min. PCR products were 
electrophoresed on 296 agarose gel, using 0.596 TBE buffer. 
Fragment sizes were determined by gel imaging analysis 
based on DNA Marker (®X174-Hinc |! digest). Finally 
genotypes were exported to Excel tables for data analysis (An 
et al, 2012). 


Data analysis 
The allele number (N4) and observed (Ho) and expected (Hz) 
heterozygosity were analyzed using Popgene 1.32 (Yeh et al, 
1999). Allele richness (Ar) was calculated using FSTAT 2.9.3 
(Hered, 1995). Since allele number is influenced by sample 
size, we used allele richness for comparison (Yan & Zhang, 
2004). Deviations from Hardy-Weinberg equilibrium (HWE) 
and linkage disequilibrium were estimated using Genepop 4.2 
(Rousset, 2008). Meanwhile the Bonferroni correction was 
conducted using SPSS 18.0. The presence of null alleles was 
detected using Micro-checker 2.2.3 (Van Oosterhout et al, 
2004). Polymorphism information content (P/C) was then 
confirmed using Microsatellite Toolkit (Zhang et al, 2010c). 
The F-statistics (Fis, Fst, and Fit) and gene flow (Nm) were 
calculated by Genetix 4.05. AMOVA was conducted using 
Arlequin 3.1 to estimate genetic variation within and between 
populations as well as for individuals (Excoffier et al, 2005). 
Popgene 1.32 was used to calculate Nei’s unbiased genetic 
distance between populations (Nei, 1978). A UPGMA system 
evolutionary tree was constructed using MEGA 5.05. 
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RESULTS 


A total of 15 polymorphic microsatellite loci were detected 
from 107 candidate primer pairs in the four tested stocks of L. 
leai via cross-species amplification. The values of Na, Ar, Ho, 
He, PIC and P for testing HWE (Px-w) at each locus in each 
stock are presented in Table 1. Totally, 114 alleles were 
detected. The allele number at each locus ranged from 2 to 
21 (mean: 7.6). Overall, specimens from HZ showed the 
highest average allele richness (7.743). PIC values for the 
four stocks ranged from 0.374 to 0.927 (mean: 0.907). 
Collectively, stocks from DT (0.412) and AQ (0.455) had the 
lowest and highest Ho, respectively, while AQ had the lowest 
He (0.791) whereas HZ had a relatively high He value (0.868). 
When the four stocks were treated as one population, no 
significant linkage disequilibrium among the loci was detected 
(P>0.05), though 12 loci showed significant (P«0.05) or highly 
significant (P«0.01) deviations from HWE (Table 1). 

F-statistics for 15 microsatellite loci among all four L. /eai 
stocks placed the mean values of Fis, Fst and Fit at 0.462, 
0.523 and 0.114, respectively. Pairwise comparisons revealed 
that the Fst ranged from 0.073 to 0.146 (0.05«Fst«0.15), 
indicating a medium differentiation among the four 
populations that was moderately closer to no differentiation 
(Wright, 1965). The number of migrats per generation (Nm) 
ranged from 1.46 to 3.18 (Table 2). 

AMOVA analysis showed that most (44.68%) of the genetic 
variation originated within individuals with only 12.5696 
variation between the four tested populations (Table 3). The 
genetic distance matrix data indicated that PY and HZ 
populations had the smallest genetic distance (0.192) and the 
highest genetic similarity (0.563). Likewise, populations from 
DT and PY showed the greatest genetic distance (0.890) and 
lowest genetic similarity (0.287) (Table 4). The UPGMA 
dendrogram further showed that PY and HZ grouped together 
and then gathered with DT and AQ, which clustered together 
(Figure 2). 


Poyang 
Hung-tse 
Dongting 
Anqing 


E O) 
0.2 0.0 


Figure 2 UPGMA dendrogram based on Nei’s genetic distances 
among Lamprotula leai stocks 


DISCUSSION 


Genetic diversity within populations 

Higher levels of genetic diversity among intensively cultivated 
stocks of aquatic species typically grants enhanced 
evolutionary potential for dealing with enviromental change, 
artificial selection, and pathogen infection (Liu & Yao, 2013; 
Wu et al, 2013), thereby improving yields and minimizing 
economic losses associated with these factors. In the present 
study, we wanted to provide basic data on the genetic 


36 WWW.ZOOres.ac.cn 


diversity and stock structure of L. leai using microsatellite 
markers to provide a foundaion for a more comprehensive 
genetic resource protection and genetic management. 
Previous studies found that sample size do not always 
correlate with expected heterozygosity, which can be a 
comparison parameter (Maudet et al, 2002; Sun, 1996). 
Compared with other mollusks, L. leai showed a higher than 
expected heterozygosity. Previously, Li et al (2009) used eight 
microsatellite markers to investigate the genetic status of H. 
cumingii, with an expected heterozygosity of H. cumingii in PY 
at 0.706, which was less than the value for L. leai (0.835) that 
we arrived at in this study. Shen et al. (2013) had previously 
isolated 14 polymorphic loci for Mytilus coruscus with an 
expected heterozygosity averaged at 0.82, while Xu et al 
(2011) estimated the mean expected heterozygosity 0.683 
and 0.759 using magnetic bead hybridization and 5'-anchored 
PCR methods to isolate microsatellite markers of L. leai. The 
disparity in these disparate findings suggests that cross- 
amplification may be more effective as compared with other 
methods. This supposition is supported by our present results, 
we showed that the overall level of heterozygosity was high 
and the PIC of 15 loci in this study showed high polymorphism 
(PIC>0.5). Together, these findings indicated that the stocks of 
L. leai have abundant genetic diversity (Botstein et al, 1980), 
but also that microsatellite loci can be used to analyze the 
genetic diversity and structure of mollusks like L. leai. More 
specifically, our findings showed that specimens from Hung- 
tse Lake had the highest Ar, average He and PIC, suggesting 
specimens from this location have a potentially larger amount 
of genetic resources, making them ideal candidate for future 
breeding efforts in China. 


Genetic differentiation among populations 

Genetic differentiation is commonly measured by Fst, Nm and 
genetic distance. Pairwise comparisons revealed medium 
genetic differentiation among the populations. AMOVA 
analysis also showed a 12.56% variation existed among these 
stocks (Wang et al, 2011b). Genetic differentiation can be 
caused by several different factors, including migration, 
genetic drift and gene mutation. Since the Nm values found in 
this study were all greater than 1, it is likely that genetic drift is 
not the main factor among L. /eai (Slatkin, 1985). However, it 
seems more likely that this phenomenon may have been 
caused by short-time farming, as intensive farming could 
affect the variation among stocks. 

The UPGMA cluster indicated that DT and AQ clustered 
together. One potential reason for this patter is that individuals 
of L. leai in Anqing Lake may have initially came from 
Dongting Lake. This possibility caused by either manual 
intervention or species migration. Study showed that genetic 
communication may related to geographic location (Ma et al, 
2007). Here Poyang Lake in the middle and lower reaches of 
the Yangtze River clustered with Hung-tse Lake, which is in 
the downstream of Huaihe River (Xu et al, 2013). These two 
lakes in close geographic proximity had similar genetic 
relationships. However, this conclusion remains to be more 
verification. 
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Table 2 
leai stocks estimated from 15 microsatellite loci 


Pairwise Fs; (below the diagonal) and the number of migrants per generation, Nm (above the diagonal) among four Lamprotula 





PY DT HZ AQ 
PY 1.46 3.18 1.47 
DT 0.146* 1.75 1.87 
HZ 0.073* 0.125* 1.52 
AQ 0.146* 0.118* 0.142* 


PY, Poyang Lake; DT, Dongting Lake; HZ, Hung-tse Lake; AQ, Anqing Lake. *: P«0.05. 


Table3 AMOVA of 15 microsatellites in the four Lamprotula leai stocks 








Source of variation df Sum of squares Variance components Percentage of variation 
Among stocks 3 187.47 0.90 12.56* 
Among individuals within stocks 115 1065.96 3.04 42.76* 
Within individuals 119 378.50 3.18 44.68* 
Total 237 1631.93 7.12 100 
*: P« 0.05. 


Table4 Neis unbiased genetic identity (above diagonal) and 
genetic distance (below diagonal) 





PY DT HZ AQ 
PY 0.287 0.563 0.361 
DT 0.890 0.398 0.574 
HZ 0.438 0.578 0.298 
AQ 0.672 0.192 0.768 


PY, Poyang Lake; DT, Dongting Lake; HZ, Hung-tse Lake; AQ, Anqing 
Lake. 


Hardy-Weinberg equilibrium 

A total of 12 loci significantly deviated from HWE after 
Bonferroni correction. Departure from HWE can be attributed 
to random genetic drift, inbreeding, and null alleles which 
contribute to heterozygote deficiency (Nei, 1987; Zheng et al, 
2009). Using Micro-checker, we found five loci including 
GQ302635, HCM08, HCM29, MP19 and APS45 with null 
alleles. This phenomenon is not unusual when analyzing 
conservation genetics using cross-amplification methods 
(Harper et al, 2003). However, each locus which departured 
from HWE can amplify at least one allele in all samples which 
means that the frequency of null alleles were not enough to 
affect the analysis (Goodman et al, 2001). In this study, Fis 
was positive and both Ho and Hg had large differences, 
suggesting an ongoing inbreeding among the tested stocks of 
L. leai. If accurate, then inbreeding is likely a critical factor 
accounting for departure from HWE, highlighting the need for 
better management and planned breeding. 


Feasibility of cross-species amplification in shellfish 

To date, 18 microsatellite loci for L. leai have been screened, 
but developing further microsatellite loci is essential for further 
more targeted research and for more cross-species 
applications. The primers derived from the relatively 
conservative flanking sequences can be used for amplification 
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across species. Some studies in shellfish have shown the 
feasibility of cross-species amplification method. For instance, 
Xu et al (2011) isolated and characterized 18 loci in L. leai, 10 
of which were successfully amplified in three H. cumingii 
populations. Wang et al (2006) confirmed that 13 loci from 32 
polymorphic microsatellite primers of Crassostrea gigas can 
amplify specific products in H. cumingii. Hai et al (2009) 
indicated that 17 primers of Perca schrenkii were amplified in 
P. fluviatilis and P. flavescens, among which 10 pairs showed 
versatility in the same genus. The microsatellite loci isolated 
in this study may then provide useful information for further 
inquiries into shellfish genetic resource information collecting 
using cross-amplification. 


CONCLUSION 


In conclusion, the four tested stocks of L. leai examined in this 
study showed both high genetic diversity and medium genetic 
differentiation, but we also encountered evidence suggesting 
ongoing inbreeding among the stocks. These results may be 
related to the life habits and reproductive characteristics of L. 
leai, or to the currently employed methods of farming and 
aquaculture. Given this result, genetic management should be 
carried out in order to maintain the genetic integrity before the 
situation deteriorates (Jia et al, 2012). Based on our results, 
we suggest that specimens from Hung-tse Lake may be 
viable parents for systematic breeding and hybrid efforts that 
may result in offspring that could offset the apparent 
inbreeding and maintain the genetic diversity we observed. 
On the whole, however, our results and methodology may be 
useful in identifying growth traits associated markers, 
constructing of genetic linkage maps, and marker-assisted 
breeding of L. leai. 
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ABSTRACT 


Multidrug resistant (MDR) pathogen infections are 
serious threats to hospitalized patients because of 
the limited therapeutic options. A novel group of 
antibiotic candidates, antimicrobial peptides (AMPs), 
have recently shown powerful activities against both 
Gram-negative and Gram-positive bacteria. 
Unfortunately, the viability of using these AMPs in 
clinical settings remains to be seen, since most still 
need to be evaluated prior to clinical trials and not all 
of AMPs are potent against MDR clinical isolates. To 
find a connection between the characteristics of 
several of these AMPs and their effects against MDR 
pathogens, we selected 14 AMPs of animal origin 
with typical structures and evaluated their in vitro 
activities against clinical strains of extensive drug- 
resistant Acinetobacter baumannii, methicillin- 
resistant Staphylococcus aureus, extended spectrum 
B-lactamase-producing Pseudomonas aeruginosa 
and extended spectrum  B-lactamase-producing 
Escherichia coli. Our results showed that these 
peptides’ hydrophilic/hydrophobic characteristics, 
rather than their secondary structures, may explain 
their antibacterial effects on these clinical isolates. 
Peptides that are amphipathic along the longitudinal 
direction seemed to be effective against Gram- 
negative pathogens, while peptides with hydrophilic 
terminals separated by a hydrophobic intermediate 
section appeared to be effective against both Gram- 
negative and Gram-positive pathogens. Among 
these, cathelicidin-BF was found to inhibit all of the 
Gram-negative pathogens tested at dosages of no 
more than 16 mg/L, killing a pandrug-resistant A. 
baumannii strain within 2 h at 4xMICs and 4 h at 
2xMICs. Tachyplesin 111 was also found capable of 
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inhibiting all Gram-negative and Gram-positive 
pathogens tested at no more than 16 mg/L, and 
similarly killed the same A. baumannii strain within 4 
h at 4xMICs and 2xMICs. These results suggest that 
both cathelicidin-BF and tachyplesin Ill are likely 
viable targets for the development of AMPs for 
clinical uses. 


Keywords: — Hydrophilic/hydrophobic ^ character; 
Multidrug resistant clinical isolate; Cathelicidin-BF; 
Tachyplesin III 


INTRODUCTION 


Multidrug resistant (MDR) pathogens continue to pose serious 
threats to hospitalized patients because there are limited 
effective therapies capable of combatting the infection. Among 
these pathogens, Acinetobacter baumannii is particularly 
intractable because several of its clinical isolates have gained 
resistance to nearly all currently available antibiotics, leading it 
to be described as "extensive drug resistant" (XDR) (Falagas & 
Karageorgopoulos, 2008). Though work on developing novel 
cocktails of antibiotics that can be used in tandem may help in 
the short-run, a long-term solution is urgently needed. 

One promising candidate source of novel antibiotic 
treatments, antimicrobial peptides (AMPs), have gained 
increased attention. These AMPs, which are innate immune 
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molecules that are widely distributed among animals, have 
previously exhibited powerful killing effects to both Gram- 
negative and Gram-positive bacteria. Unfortunately, clinical 
applications of these AMPs have made relatively little 
progress over the last decade and have been rarely reported 
(Hancock & Sahl, 2006; Lipsky et al, 2008; Vaara, 2009). The 
reasons for this lack of progress are multifaceted; aside from 
the instability of these peptides in vivo, three further reasons 
warrant some explanation. First, not all of the currently AMPs 
reported are potent AMPs. Traditionally, AMPs are purified 
directly from tissues following their activities, which may 
assure their antimicrobial activities. Currently, a growing 
number AMPs are simply predicted via bioinformatics, and 
then sometimes synthesized as peptides to test for 
antimicrobial activities. Typically, these synthetic AMPs are 
found to exhibit comparatively weak antimicrobial effects. 
Second, the structures of AMPs are highly diversified. It is 
nearly impossible to test the antimicrobial activities of these 
AMPs one by one, largely due to the intensive time and 
financial requirements. To get around this obstacle, it is 
usually necessary to select typical AMPs and then estimate 
the antibacterial activities of their corresponding groups, of 
which four predominate: helix, beta-sheet formed with 2-3 
disulfide bridges, linear peptides rich in special amino acids 
and loop peptides formed by one disulfide bridge (Vaara, 
2009). Third, the primary laboratory screening procedure of 
AMPs usually involves testing their antibacterial activities on 
several type culture strains (e.g., American Type Culture 
Collection, or ATCC strains) or several randomly selected 
clinical isolates. Though occasionally insightful, these 
findings are rarely translatable to combatting clinical MDR 
strains. 

To date, more than 2 300 AMPs have been discovered, and 
more are constantly being reported (http://aps. unmc.edu/ 
AP/main.php). This growing pool of usable materials is quite 
timely (Wang et al, 2009). In the present study, we have 
attempted and elucidate some relationship between the 
characters of several discovered AMPs and their activities 
against MDR clinical pathogens, so as to provide clues to the 
screening procedure of clinically valuable candidates. Here, 
we selected 14 AMPs of animal origin with typical structures 
to test their in vitro effects on XDR A. baumannii, methicillin- 
resistant Staphylococcus aureus (MRSA), extended spectrum 
B-lactamase (ESBL)-producing Pseudomonas aeruginosa, 
and ESBL-producing Escherichia coli. 


MATERIALS AND METHODS 


Ethics 

The protocols of this study were approved by the ethics 
committee of Kunming Medical University and ethics 
committee of Institute of Medical Biology, Chinese Academy 
of Medical Sciences and Peking Union Medical College. 
Isolation of multidrug resistant clinical strains were undertaken 
with the informed and written consent of each patient. The 
study methodologies conformed to the standards set by the 
Declaration of Helsinki and all other relevant national and 
international regulations. 
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Bacterial strains 

XDR A. baumannii strains Ab5753 and Ab5755 were isolated 
from sputum samples of two different ICU patients in the First 
Affiliated Hospital of Kunming Medical University. All other 
strains were isolated from patients at the Second Affiliated 
Hospital of Kunming Medical University. These strains tested 
in this study include the following: (I) XDR A. baumannii strain 
Ab1408 isolated from the sputum sample of one oncology 
department patient; (11) MRSA strain Sa1390 isolated from the 
sputum sample of one surgery intensive care unit (SICU) 
patient; (111) ESBL P. aeruginosa strains Pa1409 and Pa4216 
isolated from the sputum samples of two different oncology 
department patients; and (IV) ESBL E. coli strain Ec513 
isolated from the blood sample of one general surgery 
department patient. Bacterial identification was performed 
using a Vitek 32 system (bioMérieux, France). 


Antimicrobial peptides and antibiotics 

All tested AMPs were supplied by GL Biochem Ltd (Shanghai, 
China), and had a purity 295%. The tested AMPs were 
dissolved in water at 2 mg/mL and stored at -20 °C before 
use. Determining whether the AMPs are amidated on the C- 
terminus or not was made through the previously published 
literature (Table 1). Antibiotics used as references were 
supplied by Sigma-Aldrich, and dissolved in water at 2 mg/mL 
prior use, and included colistin sulfate salt (colistin), 
vancomycin hydrochloride hydrate (vancomycin), and 
tigecycline hydrate (tigecycline). 


Susceptibility testing 

The minimal inhibitory concentrations (MICs) of all AMPs and 
antibiotic references to clinical isolates were determined by 
broth dilution method in MHII following Clinical and Laboratory 
Standards Institute (CLSI) recommendations (Clinical and 
Laboratory Standards Institute, 2013). 


Time-killing curves 

Time-kill studies of tachyplesinlll, cathelicidin-BF, colistin, and 
tigecycline on XDR A. baumannii strain Ab1408 with initial 
inocula between 1x10° and 1x10’ CFU/mL were performed 
using 4xMICs, 2xMICs, 1xMIC, and 0.5xMIC concentrations. 
Samples were taken at 0, 0.5, 2, 4, 8, and 24 h after incub- 
ation. The effects of drug carryover were addressed via three 
dilution steps. Only plates with between 30 and 300 colonies 
were counted. An antibiotic was considered bactericidal when 
a reduction of 3 logio CFU/mL was achieved, as compared 
with the initial inocula (Isenberg, 2004). These tests were 
each performed in duplicate. 


Antimicrobial peptide structure prediction 

Human LL-37 and beta-defensin-2—both of which are 
longer than 30 amino acids—were sourced directly from 
NCBI PDB database PDB 2K60 and PDB 1FD3. The other 
12 AMPs were predicted by the structure prediction 
software PEP-FOLD at http://bioserv.rpbs. univ-paris- 
diderot.fr/PEP-FOLD/ (Thévenet et al, 2012), and model 1 
of each peptide with the best conformation was selected 
for further analysis. 





Table 1 Peptides tested for antimicrobial effects against drug resistant bacteria 
Name Sequence Structure Source and references 
Hyaloph 1 ilk th)*-Api: life b 
CA(1-7)M(2-9)  KWKLFKKIGAVLKVL-NH2 Helix yalophora Becropia (silt Maa Apis melon. (000 
venom) (Giacometti et al, 2003) 
[E4kjAlyteserin-1c GLKKIFKAGLGSLVKGIAAHVAS-NH2 Helix Alytes obstetricans (midwife toad) (Conlon et al, 2009) 
; Lithobates septentrionalis (mink frog) (Conlon et al, 
[DAK]B2RP GIWKTIKSMGKVFAGKILQNL-NH2 Helix 2010) 
Cathelicidin-BF —— KFFRKLKKSVKKRAKEFFKKPRVIGVSIPF Helix Bungarus fasciatus (banded krait) (Wang et al, 2008) 
LLGDFFRKSKEKIGKEFKRIVQRIKDFLRN : ; em 
LL-37 LVPRTES Helix Homo sapiens (Thomas-Virnig et al, 2009) 
GLLRKGGEKIGEKLKKIGQKIKNFFQKLVP 
Cramp QPEQ 이 S Helix Mus musculus (house mouse) (Chromek et al, 2006) 
Odo VDKPPYLPRPRPPRRIYNR-NH; ringan Proline/Arginine Oncopeltus fasciatus (milkweed bug) (Knappe et al, 
rich 2010) 
Indolicidin ILPWKWPWWPWRR-NH2 Linear, Tryptophane rich Bos taurus (cattle) (Selsted et al, 1992) 
Histatin AKRHHGYKRKFH-NH2 Linear, Histatine rich Homo sapiens (Giacometti et al, 2005) 
. " Loop, one Podisus maculiventris (spined soldier bug) (Pagès et al, 
Thanatin GSKKPVPIIYCNRRTGKCQRM M : 
disulfide bridge 2003) 
Loop, one 


Ranalexin-1Ca 


Tachyplesin III 


Beta-Defensin-2 


Alpha-defensin-2 


FLGGLMKAFPALICAVTKKC? 


KWCFRVCYRGICYRKCR-NH2* 


GIGDPVTCLKSGAICHPVFCPRRYKQIG- 
TCGLPGTKCCKKP? 
LRDLVCYCRTRGCKRRERMNGTCRKGH 


disulfide bridge 
Beta-sheet, two 
disulfide bridges 
Beta-sheet, three 
disulfide bridges 
Beta-sheet, three 
disulfide bridges 


Rana clamitans (green frog) (Halverson et al, 2000) 


Tachypleus gigas (Southeast Asian horseshoe crab) 
(Cirioni et al, 2007) 


Homo sapiens (Routsias et al, 2010) 


Mus musculus (house mouse) (Ouellette et al, 1992) 


LMYTLCCR? 





ê: For antimicrobial peptides with disulfide bridges, cysteine (C) with the same type font or underlined formed one disulfide. 


RESULTS 


Minimal inhibitory concentrations (MICs) 

All MICs are listed in Table 2. We found that all three A. 
baumannii strains were resistant to colistin according to CLSI 
standards (i.e., 24 mg/L). According to the susceptibility/ 
resistance breakpoints of tigecycline, as interpreted by the 
European Committee on Antimicrobial Susceptibility Testing 
(EUCAST) (susceptible, <1 mg/L; resistant, 24 mg/L) (The 
European Committee on Antimicrobial Susceptibility Testing 
(EUCAST) Steering Committee, 2006), Ab5753 has 
intermediate resistance to tigecycline (KDR) whereas Ab5755 
and Ab1408 are resistant to tigecycline (Pandrug resistance, 
PDR). Four AMPs, including CA(1-7)M(2-9), [D4K] B2RP, 
cathelicidin-BF, and Tachyplesin 111, showed MICs between 4 
and 16 mg/L to all of the three A. baumannii strains. Given 
that the molecular weights of these peptides (from 1 770 to 3 
638) are about three times to six times that of tigecycline 
(586), these peptides are considered as effective as 
tigecycline in inhibiting XDR A. baumannii. 

ESBL E. coli and P. aeruginosa are still sensitive to colistin, 
but the two P. aeruginosa strains were resistant to tigecycline. 
Four AMPs, namely, CA(1-7)M(2-9), [DAK]B2RP, cathelicidin- 
BF, and tachyplesinlll, showed MICs of between 4 and 8 mg/L 
to Ec513. By contrast, only two AMPs, namely, cathelicidin-BF 


and tachyplesin Ill, showed MICs between 4 and 8 mg/L to 
Pa4216 and Pa1409. This phenomenon is consistent with the 
results on tigecycline; that is, ESBL P. aeruginosa is more 
resistant than ESBL E. coli is. 

MRSA strain Sal1390 is sensitive to vancomycin, but 
intermediately resistant to tigecycline. Three AMPs, namely, 
CA(1-7)M(2-9), [D4K]B2RP, and tachyplesin 111, showed MICs 
of 16 mg/L to this MRSA strain. Notably, [E4K]Alyteserin-1c 
and cathelicidin-BF, which showed a generally effective 
inhibitory activity on all MDR Gram-negative bacteria, had no 
effect on Gram-positive MRSA. 


Time-killing kinetics 

Figure 1 shows the time-killing kinetics. Cathelicidin-BF kills 
XDR A. baumannii strain Ab1408 within 4 h at 2xMICs and 2 
h at 4xMICs respectively. This is similar to the action of 
colistin, i.e., Cathelicidin-BF was also able to rapidly kill 
bacteria. By comparison, tachyplesin Ill kills XDR A. 
baumannii after 4 h of incubation at 2xMICs and 4xMICs. 
Compared with colistin, tigecycline showed a markedly 
slower killing effect (after 8 h of incubation at 2xMICs and 
4xMICs), but this effect lasted longer (24 h at 4xMICs). 
The optimum killing effects of cathelicidin-BF and tachyp- 
16910 111 appeared to be 4 h of incubation, and afterward 
gradually declined. 
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Time-killing curves of colistin 
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Figure 1 Time-killing curves of colistin, tigecycline, cathelicidin-BF, and tachyplesin 111 on one XDR Acinetobacter baumannii strain 


(Ab1408) 


Initial inocula between 1x10? and 1x10’ CFU/mL were performed using 4xMICs, 2xMICs, 1xMIC, and 0.5xMIC concentrations. Samples were obtained at 0, 


0.5, 2, 4, 8, and 24 h after incubation. 


Amphipathic structure distributions of AMPs 

Although CA(1-7)M(2-9), [DAK]B2RP, [E4K]Alyteserin-1c, 
cathelicidin-BF, LL-37, and cramp are all classified as helical 
peptides by their secondary structures, these AMPs showed 
different antibacterial effects. Since the amphipathic structure 
distributions are generally considered crucial for AMPs to kill 
bacteria (Brogden, 2005), we studied the three-dimensional 
structures of the eight AMPs that were found to be effective 
on our tested MDR clinical isolates (Figure 2). [E4K]Alyte- 
serin-1c and cathelicidin-BF, which were effective on all three 
Gram-negative bacteria, both have classical long linear 
amphipathic structures, withhydrophobic regions on one side 
and hydrophilic regions on the other side along the linear 
peptides. Similar structures exist in LL-37 and cramp, but 
these AMPs only showed slight activity on a portion of these 
Gram-negative strains. The helical peptides of CA(1-7)M(2-9) 
and [D4K]B2RP, which could inhibit not only Gram-negative 
but also Gram-positive bacteria, tend to manifest a more 
contractive style, with hydrophilic terminals separated by a 
hydrophobic intermediate section. A similar amphipathic 
structure distribution can be seen in tachyplesin Ill, which is 
also effective on both Gram-negative and Gram-positive 
clinical isolates, but is classified as a beta-sheet formed by 
two disulfide bridges according to its secondary structure. 


44 WWW.ZOOres.ac.cn 


Interestingly, a similar structure also exists in Ranalexin-1Ca, 
which only showed slight activity to Gram-positive strains, but 
was ineffective on Gram-negative strains. 

AMPs (mainly beta-sheet AMPs like beta-defensin-2 and 
alpha-defensin-2 that are formed with 3 disulfide bridges, 
and linear AMPs like oncocin, indolicidin, histatin and 
thanatin, which are rich in special amino acids) with neither 
classical long linear amphipathic structures nor hydrophilic 
terminals separated by hydrophobic intermediate sections, 
showed no effects on the tested MDR clinical isolates (Table 
2, Figure 3). 


DISCUSSION 


In total, more than 2 300 antimicrobial peptides have been 
reported to date, but as we mentioned earlier, clinical 
applications of these AMPs are held up by a number of 
roadblocks. In this study, we collected a set of AMPs with 
different structures to evaluate their in vitro effects on typical 
MDR clinical isolates, which included XDR A. baumannii, 
MRSA, ESBL-producing P. aeruginosa and E. coli. Our results 
showed that helical AMPs appear to be more effective at 
killing MDR bacteria than other types of AMPs. By contrast, 
linear AMPs rich in special amino acids and AMPs with 







Figure 2 Structures and three-dimensional hydrophilic 
Ihydrophobic arrangements of effective antimicrobial peptides 
The backbones of these AMPs are shown in flat ribbons, and the three- 
dimensional surfaces of these AMPs are shown in solid. Blue: hydrophilic a 
regions; Red: hydrophobic regions. A: CA(1-7)M(2-9); B: [D4K]B2RP; C: Figure 3 Structures and three-dimensional hydrophilic/ hydro- 


Tachyplesin 111; D: Ranalexin-1Ca; E: cathelicidin-BF; F: [E4K]Alyteserin-1c; phobic arrangements of ineffective antimicrobial peptides 

G: LL-37; H: cramp. Dotted lines indicate different hydrophilic/ hydrophobic The backbones of these AMPs are shown in flat ribbons, and the three- 
region arrangement styles between A-D (both Gram-negative and Gram- dimensional surfaces of these AMPs are shown in solid. Blue: hydrophilic 
positive effective, except for D, i.e., Ranalexin-1Ca, which is only Gram- regions; Red: hydrophobic regions. A: Indolicidin; B: Histatin; C: Thanatin; 
positive effective) and E-H (Gram-negative effective). D: Oncocin; E: beta-Defensin-2; F: alpha-defensin-2. 


Table2 MICs (mg/L) of AMPs and antibiotics against drug resistant bacteria 








KERAN XDR A.baumannii ESBLs E.coli & P. aeruginosa MRSA 

AMPs & antibiotics 
Ab5753 Ab5755 Ab1408 Ec513 Pa1409 Pa4216 Sa1390 

CA(1-7)M(2-9) 8 4 8 8 64 32 16 
[DAK]B2RP 8 4 16 8 32 32 16 
[E4K]Alyteserin-1c 64 8 64 32 2128 64 2128 
Cathelicidin-BF 16 4 16 8 8 4 2128 
LL-37 2128 32 2128 128 2128 2128 2128 
Cramp 2128 128 2128 128 2128 2128 2128 
Oncocin 2128 2128 2128 2128 2128 2128 2128 
Indolicidin 2128 2128 2128 2128 2128 2128 2128 
Histatin 2128 2128 2128 2128 2128 2128 2128 
Thanatin 2128 2128 2128 2128 2128 2128 2128 
Ranalexin-1Ca 2128 2128 2128 2128 2128 2128 64 
Tachyplesin III 8 8 16 4 8 8 16 
Beta-Defensin-2 2128 2128 2128 2128 2128 2128 2128 
Alpha-defensin-2 2128 2128 2128 2128 2128 2128 2128 
Colistin 2128 16 8 «1 «1 «1 2128 
Vancomycin 2128 2128 2128 2128 2128 2128 2 
Tigecycline 2 4 4 2 16 16 2 


Abbreviations: AMP: antimicrobial peptides; MIC: minimum inhibitory concentration; XDR: extensive drug resistance; ESBLs: extended spectrum beta- 
lactamase producing; MRSA: methicillin-resistant S. aureus. 
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disulfide bridges (with the exception of tachyplesin 111) are 
generally less effective in combatting MDR bacteria as 
compared with helical AMPs. Three-dimensional analysis 
showed that the hydrophilic/hydrophobic arrangements of 
these peptides, and not their secondary structures, seem to 
contribute more to their efficacy and position on the 
antimicrobial spectra. Peptides with long linear amphipathic 
structures were found to be effective only on Gram-negative 
pathogens, whereas peptides with more contractive styles 
with hydrophilic terminals separated by a hydrophobic 
intermediate section appeared to be effective on both Gram- 
negative and Gram-positive pathogens. Even though all AMPs 
are reported as antimicrobial peptides, only 2 of the 14 tested 
AMPs (tachyplesin IIl and cathelicidin-BF) were effective 
against all the tested MDR pathogens within their 
antimicrobial spectrum and with low MICs (<16 mg/L). Since 
the net charge these AMPs carried does not correlate with 
their antibacterial activities (data not shown), the underlying 
mechanisms that lead to different MICs to kill bacteria 
remains to be further determined in more targeted follow-up 
studies. 

Our results also showed that Cathelicidin-BF kills bacteria 
as quickly as colistin (within 4 h at 2xMICs and 2 h at 
4xMICs), while Tachyplesin Ill is slower (after 4 h of 
incubation at 2xMICs and 4xMICs) than colistin, but faster 
than tigecycline. Considering the difference of the amphipathic 
Structure distributions between these two peptides, the 
different rate of killing bacteria may imply a difference in the 
mechanism underlying their activities. Notably, the killing 
effects of these two potent peptides all declined after 4 h of 
incubation, which may, in part, be due to the nature of the 
peptides. Again, further targeted studies are needed to shed 
more light into the differences of these AMPs. 

In conclusion, in the present study not all of the reported 
AMPs were effective against several tested MDR clinical 
isolates. Of the 14 potential AMPs, only two, Tachyplesin 111 
and cathelicidin-BF, which differ in both their secondary 
structures and three-dimensional hydrophilic/ hydrophobic 
arrangements, showed potent activities (<16 mg/L) against 
nearly all of the MDR clinical isolates within corresponding 
antimicrobial spectrum in vitro. What differentiates these two 
AMPs in terms of their efficacy to kill bacteria compared with 
other AMPs with similar hydrophilic/ hydrophobic arran- 
gements remains to be determined. 
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ABSTRACT 


To explore the neural mechanisms mediating aging- 
related visual function declines, we compared the 
expressions of brain-derived neurotrophic factor 
(BDNF) and its high affinity receptor-tyrosine kinase 
B (TrkB) between young and old adult cats. Nissl 
staining was used to display neurons in each layer of 
the lateral geniculate nucleus (LGN). The BDNF- and 
TrkB receptor-immunoreactive neurons were labeled 
immunohistochemically, observed under optical 
microscope and photographed. Their neuronal 
density and immunoreactive intensity were measured. 
Results showed that the mean density of the Nissl 
stained neurons in each LGN layer were comparable 
between old and young adult cats, and their BDNF 
and TrkB proteins were widely expressed in all LGN 
layers. However, compared with young adult cats, 
both the density and optical absorbance intensity of 
BDNF- and TrkB-immunoreactive cells in each LGN 
layer in old cats were significantly decreased. These 
findings indicate that the decreased expressions of 
BDNF and TrkB proteins in the LGN may be an 
important factor inducing the compromised inhibition 
in the central visual nucleus and the functional visual 
decline in senescent individuals. 


Keywords: Cat; Lateral geniculate nucleus; BDNF; 
TrkB; Age-related change 


INTRODUCTION 


Progressive recession in visual capacity, which includes 
decreased visual acuity, insensitiveness to visual contrast, 
inability to determine orientation or direction of objects and 
slowdown in processing visual information, can be found 
during the normal course of aging (Hua et al, 2006; 
Schmolesky et al, 2000). However, the neuronal mechanisms 
underlying age-related functional visual decline are still 
unclear (Hua et al, 2011; Zhou et al, 2011). How changes in 
neuronal plasticity affect age-related functional visual decline 
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and the correlated molecular or cellular mechanisms are quite 
controversial (Hua et al, 2006; Liang et al, 2012b; Schmolesky 
et al, 2000; Zhou et al, 2013). For example, previous electr- 
ophysiological studies reported that the reactive characters of 
neurons in the lateral geniculate nucleus (LGN) were not 
significantly correlated with age and that aging might affect 
neuronal functions of the visual cortex, especial the advanced 
visual cortex (Hua et al, 2006; Liang et al, 2012a, 2012b; 
Zhou et al, 2011). However, these results could not rule out 
the possibility that both the morphology and neurotransmitter 
systems had changed during aging (Vidal et al, 2004), and 
these changes could influence the reactive characters of the 
visual cortical neurons to visual stimuli through synaptic 
connections. Other studies reported that aging was also 
accompanied with decreased expression of inhibitory 
neurotransmitters (Such as GABA) in the visual cortex (Hua 
et al, 2008), which might induce the insensitiveness and 
decline in selectivity of the visual cortical neurons to visual 
stimuli (Hua et al, 2006; Leventhal et al, 2003; Zhou et al, 
2011). Studies on neurotrophic factors found that brain- 
derived neurotrophic factors (BDNF) are vital in the 
development, survival, and lesion repair of neurons and 
neuronal plasticity via receptor mediated signal transductions 
(Ichim et al, 2012; Kwon et al, 2011; Numakawa et al, 2010). 
Some reports have shown that the BDNF signal system can 
modulate activity-dependent intracortical inhibitory (Rutherford 
et al, 1997) and promote the synthesis and transportation of 
GABA (Sánchez-Huertas & Rico, 2011; Vaz et al, 2011; 
Waterhouse et al, 2012). However, although a few studies 
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reported that age-related changes were found in the cerebral 
expressions of BDNF and receptors (Erickson et al, 2010; 
Hayakawa et al, 2008; Luellen et al, 2007; Silhol et al, 2005), 
it is too early to conclude that the decreased expression of 
BDNF and its receptors during senescence is the reason for a 
decline of inhibitory neurotransmitter GABA. 

Our previous immunohistochemistry results showed that the 
expression of GABA in the dorsal LGN (dLGN) decreased 
significantly during aging (Tong et al, 2006); however, whether 
this phenomena is correlated with the expression down- 
regulation of BDNF and its receptors remains unclear. To 
understand the cellular and molecular mechanisms of age-related 
functional visual decline, we explored the expression differences 
of BDNF and its high affinitive receptor, tyrosine kinase receptor 
B (TrkB), in the dLGN of young and old adult cats. 


MATERIALS AND METHODS 


Animal subjects and experimental reagents 

Four young (1-3 years old, 2-3.5 kg) and four old (10-13 years 
old, 2-3.5 kg) adult cats were examined ophthalmoscopically 
prior to experimentation to confirm that no optical or retinal 
problems impaired their visual function. All experiment 
procedures were performed strictly in accordance with the 
guidelines published in the NIH Guide for the Care and Use of 
Laboratory Animals. After a deep anesthesia by ketamine HCI 
(40 mg/kg, im) injection, cardiac perfusion was immediately 
performed on the animals (0.9% saline water). When pale 
livers were observed, a pre-fixture with 0.1 mol/L PBS (10% 
formalin, 2.5% glutaraldehyde, pH 7.2-7.4, 200 mL/kg) was 
performed. The isolated brainstem was fixed for 2 h and the 
LGN was then separated from it and transferred into a fixing 
solution (10% formalin, 2.5% glutaraldehyde, 30% sucrose) till 
the tissue sank to the bottom. Every four continuous coronal 
frozen slices (30 um) were selected as a group for Nissl 
staining, BDNF and TrkB immunohistochemistry labels and 
negative controls. The slice groups were separated from each 
other by an interval of five slices. Every cat included 10 
groups of slices. 


Neuronal Nissl staining 

Slices were stained in 0.1% cresyl violet solution for 5 min at 
room temperature, rinsed with distilled water, dehydrated by 
gradient alcohol, transparentized by xylene and then sealed 
with gum. Stained slices were used to determine the layering 
structures of the LGN and the neuronal density of each layer. 


Immunohistochemical labeling 

The frozen slices were incubated in 3% H202 for 15 min at 
room temperature to eliminate endogenous peroxidase 
activity and were then rinsed with distilled water, incubated in 
0.3% TritonX-100 PBS for 20 min at room temperature, 
incubated in solution with 5% fetal bovine serum protein for 20 
min at room temperature to seal the specific reactive areas, 
incubated in rabbit-anti-mouse BDNF and rabbit-anti-human 
TrkB multiclonal antibodies (1:100, primary antibodies) for 36 
hours at 4°C and then rinsed three times with PBS for 5 min, 


incubated in biotinylated secondary antibody (goat-ant-rabbit) 
for 20 min at room temperature and then rinsed three times 
with PBS for 5 min, incubated in streptavidin-biotin complex 
(SABC) working solution for 20 min at room temperature and 
then rinsed three times with PBS for 5 min, developed by DAB 
to produce colorimetric end products, dehydrated by gradient 
(80%, 95% and 100%) alcohol, transparentized by xylene and 
then sealed with neutral gum. The only difference in the 
treatment of the negative control slices was replacing 
primary antibodies with PBS. All immunohistochemistry kits 
and the DAB substrate were products of Boster, Wuhan, 
China. 


Statistical analysis 

Nissl stained slices and BDNF and TrkB immuno-reactive 
slices were observed under a microscope (Olympus BX-51). 
Images were collected by Image-Pro Express 6.0 software. All 
related morphological parameters went through quantitative 
analysis, including the Nissl stained neuronal densities in 
each LGN layer, densities of BDNF and TrkB immuno-reactive 
positive neurons and intensities of immuno-reactions 
(evaluated by average optical absorbance, where higher 
absorbance values indicated stronger immuno-reactions). 
Data analyses were performed unbiased. 

Images of slices were initially acquired under 40x ocular 
magnification and the layering structures (layer A, layer A; and 
layer C) of the LGN were then discriminated under 100x 
ocular magnification. Ten view fields (50 umx50 um) of each 
layer were randomly selected under 400x ocular magnification 
to calculate the densities of Nissl stained neurons and BDNF 
and TrkB immuno-positive neurons (cells/mm?), respectively. 
A typical Nissl stained neuron was characterized with royal 
purple Nissl bodies and a clear nucleus in soma, and a typical 
BDNF/TrkB immuno-positive neuron was characterized with 
obvious immuno-positive matter and a clear nucleus in soma. 
The optical absorbance values of 20 immuno-positive neurons 
from each slice were randomly selected by Image-Pro 
Express 6.0 software, and the mean value was taken as the 
index indicating the intensity of immuno-reactions. All data 
were analyzed via SPSS 13.0 software and expressed as 
means+SD, with P<0.05 being considered statistically 
significant. 


RESULTS 


Densities of Nissl stained neurons 

A three-layer-structure of the LGN (layer A, layer A; and layer 
C) (Figure 1) was shown in the Nissl stained slices. Blue or 
dark blue somas in various shapes, such as oval, tapered or 
polygon, were visualized in each layer. The results showed 
that age-independently, the neuronal density of layer A was 
significantly higher than that of layer A4 (Fa,158)=31.04, 
P«0.001), and neuronal densities of layer A and A; were both 
higher than that of layer C (Fa15s=212.04, P<0.001; 
F1.158)=85.72, P«0.001). No differences were found in the 
neuronal densities of each layer between young and old adult 
cats (F1,238)=1.88, P=0.17, Table 1). 
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Figure 1 Nissl stained laminar organization (left, layer A, A; and C) and neurons (right) in the LGN of cats (Scale bar=50 um) 


Densities of BDNF and TrkB immuno-positive neurons 

Immunohistochemical staining results showed that BDNF and 
TrkB immuno-positive neurons with brown or dark brown 
somas and their fiber distributions were found in every layer of 
the LGN in young and old adult cats (Figure 2). Age- 
independently, the densities of BDNF immuno-positive 
neurons in layer A and A, were both significantly higher than 
that of layer C (F(1,158y727.12, P«0.001; F(1458j740.06, P«0.001, 
Table 1), whereas, those between layer A and A4 were 





Figure 2 BDNF immuno-positive neurons (A, B, C, D) and TrkB 
immuno-positive neurons (E, F, G, H) in layer A of the LGN in 
young (A, C, E, G) and old (B, D, F, H) adult cats (Scale bar=50 um) 
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comparable (F(1,158=1.33, P=0.25). However, although no 
differences were found in the densities of BDNF immuno- 
positive neurons in each layer among cats at similar ages 
(F(323ej71.11, P=0.35), those in old adult cats were 
significantly lower than those in young adult cats (Fa,238)= 
125.39, P«0.001, Table 1). 

Age-independently, the densities of TrkB immuno-positive 
neurons in layer A and A; were both significantly higher than 
that of layer C (F(1,158)=28.16, P<0.001; F,158)=29.92, P<0.001, 
Table 1), whereas, those between layer A and A, were 
comparable (F(1,15s=0.02, P=0.89). However, although no 
differences were found in the densities of BDNF immuno- 
positive neurons in each layer among cats at similar ages 
(F(3236j70.22, P=0.88), those in old adult cats were 
significantly lower than those in young adult cats (Fa,238)= 
208.54, P«0.001, Table 1). 


Immuno-intensities of BDNF and TrkB immuno-positive 
neurons 

The mean optical absorbance values of the BDNF immuno- 
positive neurons in the LGN of old cats were significantly 
lower than those in young cats (F(78;718.74, P«0.001), 
specifically, the mean absorbance values of layer A, A; and C 
were decreased by 26.4%, 20.3% and 36.3%, respectively. 
The mean absorbance values of the TrkB immuno-positive 
neurons in the LGN of old cats were also significantly lower 
than those in young cats (F(178j728.46. P«0.001), specifically, 
the mean optical absorbance values of layer A, A: and C were 
decreased by 25.696, 26.596 and 35.396, respectively (Figure 3). 
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Figure 3 Mean optical absorbance values of BDNF and TrkB 
immuno-positive neurons in the LGN of young and old adult cats 
(*: P«0.01) 


DISCUSSION 


Vision is the most important way for higher animals, including 


Table 1 


Densities of BDNF, TrkB immuno-positive neurons in the dLGN of young and old cats 


Lamina organization of LGN 





Animal subjects 


A A1 





Densities of Nissl stained neurons (cells/mm?) 





YC1 3 510.0:587.8 
YC2 2 510.0:481.8 
YC3 2 550.0+259.3 
YC4 2 850.0+406.2 
Mean 2 855.04592.7* 
001 2 700.0+434.6 
002 2 810.0+481.8 
003 2 930.0+490.0 
004 2 690.0+486.4 
Mean 2 782.54465.7* 


2 650.0+462.5 
2 430.0:432.2 
2 440.0+320.4 
2 510.0+587.7 
2 507.54451.4* 
2 380.0+482.6 
2 430.0+503.4 
2 390.0+296.1 
2 020.0+297.4 
2 305.04425.4* 


2 030.0+336.8 
1 590.0+159.5 
1 880.0+312.0 
1 920.0+234.8 
1 855.0£307.1* 
1 980.0+278.1 
1 680.0+193.2 
1 750.0+222.4 
1 830.0+231.2 
1 810.0+250.9* 





Densities of BDNF immuno-positive neurons (cells/mm?) 











YC1 2 320.0+315.6 2 640.0+430.0 2 040.0+350.2 
YC2 2 400.0+461.9 2 280.0+423.7 2 000.0+326.6 
YC3 2 320.0+253.0 2 120.0+379.5 2 160.0+337.3 
YC4 2 360.0+397.8 2 440.0+461.9 1 840.04337.3 
Mean 2 350.0+353.0** 2 370.0+451.1** 2 010.0+344.8** 
001 1 640.0£295.1 2 000.0+377.1 1 280.0£253.0 
OC2 1 760.0£337.3 1 840.04386.4 1 520.0£367.6 
OC3 1 960.04350.2 1 800.0£282.8 1 080.0+270.0 
004 1 560.04295.2 1 920.04367.6 1 240.0£227.1 
Mean 1 730.0£343.6** 1 890.000+350.7** 1 280.0£316.4** 
Density of TrkB immuno-positive neurons (cells/mm?) 
YC1 2 280.0+369.5 2 520.0+329.3 2 000.0+326.6 
YC2 2 360.0+397.8 2 280.0+388.7 2 040.0+266.7 
YC3 2 320.0+168.7 2 160.0+386.4 2 120.04423.7 
YC4 2 400.0+533.3 2 480.0+454.1 1 880.04379.5 
Mean 2 340.04379.5** 2 360.0+410.6** 2 010.0+350.8** 
001 1 600.0£266.7 1 760.0£279.7 1 200.0+188.6 
OC2 1 800.04339.9 1 560.0£227.1 1 600.0+233.0 
OC3 1 920.0£367.6 1 670.04259.7 1 120.0£315.5 
OCA 1 640.0£350.2 1 960.04+380.2 1 280.04253.0 
Mean 1 740.0£345.5** 1 737.5£311.2** 1 300.0+289.3** 


YC1-4 and OC1-4 represent four different young and four different old adult cats, respectively. *-P»0.05; **: P«0.001. 


humans, to acquire environmental information. Age-related 
functional visual decline can therefore directly lower life quality. 
Understanding neuronal networks and the molecular and 
cellular mechanisms underlying age-related functional visual 
decline is helpful in delaying the effects of senescence on 
vision and stimulating visual ability in seniors (Zhang et al, 
20083). 

In higher animals and humans, the processing of visual in- 


formation includes three different levels, that is, the retina, 
LGN and visual cortex from low to high. Previous electro- 
physiological studies indicate that aging may change the 
reactive characters of visional neurons to visional stimuli, e.g., 
the ability to determine the locations of visional stimuli and 
confirm moving directions, sensitivity to contrasts and 
selectivity in dimensional and temporal frequencies (Fu et al, 
2013; Hua et al, 2006; Leventhal et al, 2003; Liang et al, 2010; 
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Schmolesky et al, 2000; Wang et al, 2005; Yang et al, 2009; 
Zhang et al, 2008b; Zhou et al, 2011). However, although 
these studies did not show significant influences of senes- 
cence on the characteristics of neurons on the LGN (Spear et 
al, 1994; Zhou et al, 2013), it is still possible that age-related 
changes are happening at the molecular or cellular level in the 
LGN. Our earlier work found that GABA expressions in the 
LGN of old cats were significantly lower than those in young 
adult cats (Tong et al, 2006). We assume that the functional 
effects of decreasing GABA may be compensated by down- 
stream pathways or other factors, and therefore senescence 
may influence the function of LGN at the cellular level. 

Tong et al (2006) reported that during the course of aging, 
the neuronal densities of the LGN maintained relatively stable, 
while the expression of GABA was significantly decreased; 
however, the underlying mechanisms for this remain unclear. 
Studies on neurotrophic factors showed that BDNF and its 
receptor (including high affinitive Trk receptor and low high 
affinitive p75 receptor) introduced signal transductions play 
vital roles in the development, survival and synaptic plasticity 
of neurons (Ichim et al, 2012; Kwon et al, 2011; Numakawa et 
al, 2010; Wong et al, 2013), including GABAnergic neurons 
(Cheng & Yeh, 2005; Jiao et al, 2011; Rutherford et al, 1997; 
Sakata et al, 2009; Waterhouse et al, 2012), and may be 
directly involved in the synthesis and transportations of GABA 
(Sanchez-Huertas & Rico, 2011; Vaz et al, 2011). However, 
whether the age-related GABA decrease is also correlated 
with the abnormal BDNF signal transduction still needs further 
evidence. 

In the present study, no differences were found in the 
neuronal densities of each LGN layer, indicating relative 
stable neuronal quantities during the aging course, which is 
consistent with previous research (Tong et al, 2006). However, 
both the densities and intensities of the immuno-activity of 
BDNF and TrkB immuno-positive neurons in each LGN layer 
of old adult cats were significantly lower than those in young 
adult cats, which is consistent with results observed in old rats 
and the hippocampus of patients with Alzheimer's disease 
(Chapman et al, 2012; Erickson et al, 2010; Luellen et al, 
2007; Tapia-Arancibia et al, 2008), indicating that the 
expressions of BDNF and TrkB decrease and the BDNF-TrkB 
signal transductions decline during the course of aging. 

In summary, the present study indicated that senescence 
could induce the decreased expression of BDNF and TrkB 
receptors in each LGN layer, and the BDNF-TrkB signal 
transduction decline could decrease the expression of GABA. 
The weakening inhibitory effects of nerve endings on synaptic 
connections could influence the sensitivity and selectivity of 
visual neurons to visual stimuli and then induce age-related 
functional visual decline. 
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ABSTRACT 


It is well known that excessive long-term alcohol 
consumption is harmful, especially in pregnant 
women. In the present study, the Kunming white 
mouse was used as an animal model and indirect 
immunofluorescence was performed to analyze the 
toxic effects of alcohol on early pre-implantation 
embryos. H3K9 acetylation immunofluorescence 
could not be detected in MII oocytes. H3K9 
acetylation levels in the treatment group were higher 
than in the control group during the morula stage, 
and contrary to results during the blastocyst stage. 
Other stages showed no obvious differences for in 
vivo embryos. For in vitro embryos, almost no 
difference was found between the two experimental 
groups across all stages, and both groups showed 
increasing H3K9 acetylation levels (except at the 2- 
cell stage). This study shows that H3K9 acetylation 
levels in early pre-implantation embryos are notably 
impacted by excessive alcohol ingestion by females. 
These data are the first step in understanding the 
epigenetic mechanism of alcohol toxicity in early pre- 
implantation mouse embryos. 


Keywords: Alcohol intake; Pre-implantation embryos; 
H3K9 acetylation; Mice; Embryo development 


INTRODUCTION 


Long-term excessive drinking is severely harmful, especially in 
pregnant women. Fetal alcohol spectrum disorder (FASD), with 
syndromes including pre- and postnatal developmental 
retardation, craniofacial paramorphia, central nervous 
abnormalities, congenital heart diseases and mental and 
behavior disorders, can be induced by excessive drinking 
during progestation or gestational periods (Mattson et al, 2011). 
Studies into the toxicology of alcohol have mainly focused on 
the effects of excessive drinking on tissues, such as alcohol 
induced fatty liver (Li, 2003), cerebral injury (Zhao & Feng, 2012) 
and retinopathy (Jiang et al, 2007), or neonatal risks and 
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syndromes following maternal alcohol exposure (Mattson et al, 
2011). Few studies have been done on the toxicology of alcohol 
on early pre-implantation embryos in mice. Huang et al 
(2013a,b) found that alcoholism reduces germocyte activity and 
induces abnormal DNA methylation in early pre-implantation 
embryos. The aim of the present study is to understand the 
effects of alcoholism on H3K9 acetylation levels of pre- 
implantation embryos in mice from an epigenetic perspective. 
The influence of epigenetic modification on biological 
phenotype has become a research hotspot. Most congenital 
diseases are the result of genetic and environmental factors 
and genetic mutation itself may not fully explain their 
nosogenesis. Epigenetic modifications include DNA 
methylation, genomic imprinting, histone modification and 
non-encoding RNA, etc. Histones are highly alkaline proteins 
found within chromosomes that package and order DNA into 
structural units. Histone H3 and H4 feature a N-terminal tail 
protruding from the globular nucleosome core, which can 
undergo several different types of  post-translational 
modification that influence cellular processes, including 
reversible covalent modifications, such as histone acetylation, 
etc. Histone acetylation plays a vital role in the development 
of mammals and its abnormality results in chromatin 
reorganizing, and genetic transcriptional disorders in the cell 
cycle, cell differentiation and apoptosis (Vigushin & Coombes, 
2002). Histone acetyltransferase (HATs) and histone 
deacetylases (HDACs) are a pair of enzymes that regulate 
genetic acetylation and deacetylation. By adjusting the 
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electrostatic attraction and sterical hinderence between DNA 
and histone, thereafter influencing the tightness of chromos- 
omes, HATs and HDACs can repress or initiate gene 
transcription (Haycock, 2009). Histone H3K9 acetylation is 
one of the signs of gene transcription initiation (Morinobu et al, 
2004) and its abnormal acetylation can induce abnormal 
chromatin structure and abnormal gene expression. Fetal 
malformations or abnormal developments are usually the 
result of protein misexpression. For example, imbalanced 
histone acetylation can cause the misexpression of cardi- 
ogenesis-related genes and bring on obstacles during the 
differentiation of mesenchymal stem cells (MSCs) to 
myocardial cells (Lin, 1981). Fisher et al (1986) reported that 
alcohol elevates acetylation levels of histone H3K9 in 
myocardial progenitor cells and then induces misexpression 
of cardiogenesis-related genes. 

Earlier studies from our group suggest that epigenetic 
modifications of pre-implantation embryos in mice are 
extremely sensitive to environmental change. Alcohol-induced 
genomic DNA methylation in germocytes and zygotes, H3K9 
methylation and H3K27 acetylation could interfere with the 
orderly establishment of epigenetic patterns in pre- 
implantation embryos in mice (Huang et al, 2013a, 2013b). To 
better understand epigenetic patterns under the influence of 
alcohol we force-fed Kunming (KM) mice with alcohol to 
stimulate external inducing conditions of FASD, and then 
monitored the acetylation level of H3K9. 


MATERIALS AND METHODS 


Experimental animals 

Female KM mice (4—5 weeks old, 30+2 0) were provided by 
the Experimental Animal Center of Anhui Medical University 
and were housed at 24+2 °C with lights on from 0800h-2200 
h. Food and water were available ad libitum. Animals were 
allowed to adapt to the housing condition for one week. All 
animal procedures were carried out in accordance with 
Fuyang Teachers College Animal Care and Use Committee. 


Experimental reagents and equipments 

Pregnant mare serum hormone (PMSG) and human chorionic 
gonadotrophin (hCG) were products of the Second Hormone 
Factory of Ningbo, China. Saline water was from Shanghai 
Huayuan Anhui Jinhui Pharmaceutical, China. The antibody of 
H3K9 acetylation was from Epigentek Group Inc. (New York, 
USA). Secondary goat-anti-rabbit antibody, IgG-FITC, was 
from Bioss Biotechnology (Beijing, China). Petri dishes were 
from Corning and all other reagents were products of Sigma. 
Equipment included a stereomicroscope (SMZ1500, Nikon, 
Japan), CO» incubator (2323-2, Thermo, USA), electronic 
balance (AR224CN, Ohaus, USA) and laser scanning cofocal 
microscope (TCS-SP5, Leica, Germany). 


Animal model establishment and animal grouping 

Animal models were established according to previous 
publications (Dole & Gentry, 1984; Dole et al, 1985; Griffin et 
al, 2009; Hwa et al, 2011) indicating that the behavioral and 


physiological features in mice force-fed with 2096 (4—5 g/kg) 
alcohol bear a resemblance to those in human alcoholics. 
Mice were randomly separated into in vivo and in vitro groups. 
In the in vivo group, animals were further divided into an 
alcohol group and control group. The alcohol treatment 
animals (n=20) were force-fed 20% (4-5 g/kg) alcohol for 30 
days and then force-fed saline water (4-5 g/kg) for another 4 
days, ad libitum; control animals (n=20) were force-fed saline 
water 34 days. In the in vitro group, animals were divided into 
an alcohol group and control group. The alcohol treatment 
animals (n=10) were force-fed 20% (4-5 g/kg) alcohol for 30 
days ad libitum; the control animals (n=10) were force-fed 
saline water 30 days. Zygotes and embryos at different stages 
were cultured in Petri dishes for 4 days. 


Collection of MII stage ovum and pre-implantation 
embryos from in vivo and in vitro groups 

Animals were given 10 U PMSG (i.p.) at 1730h on the last day 
of force-feeding and given 10 U hCG (i.p.) 48 hours later.Five 
females were randomly selected to not mate with males and 
14 hours after hCG injection, their ampulla were punctured to 
collect MII stage ovum. The other females were mated with 
adult males (1:1). At 0800h on the day after mating, vaginal 
suppositories were checked and if the vaginal suppository 
was identified then female fertilization was confirmed. /n vivo 
flushing method (using syringe, HCZB fluid and observation 
under a stereomicroscope) was applied in the in vivo group to 
collect zygotes (at 14 hours after fertilization and zygotes can 
be obtained by puncturing ampulla directly), 2-cell-, 4-cell-, 8- 
cell-embryo, morula and blastula (at 30, 48, 60, 72 and 88 
hours after fertilization) at different stages of embryonic develo- 
pment. For the in vitro group, zygotes collected 14 hours after 
fertilization were rinsed with CZB sugar-free culture fluid and 
then cultured with CZB fluid (sugar-free) in the incubator at 37 
C, 5% CO». Forty-eight hours later, 4-cell-embryos were 
transferred into CZB fluid (with sugar) and cultured in the 
incubator at 37 'C, 5% CO». Zygotes, 2-cell-, 4-cell-, 8-cell- 
embryo, morula and blastula collected in vitro were rinsed with 
PBS (0.1% PVA), fixed with 4% paraform and saved. 


Indirect immunofluorescence detection of H3K9 acetylation 
in pre-implantation embryos at different developmental 
stages 

Fixed zygotes and embryos at different stages were rinsed 
with PBS (0.1% PVA) and then incubated in PBS with 0.2% 
TritonX-100 at 37.5 'C for 1 hour; rinsed with PBS (0.1% BSA); 
sealed with PBS (0.1% BSA) for 1-2 hours at 37.5 C; 
incubated in H3K9 acetylation antibody (1:100, diluted with 
PBS (196 BSA)) over night at 37.5 °C; rinsed for 5 min, 3 times 
with PBS (0.1% Tween 20 and 0.01% TritonX-100); incubated 
in secondary goat-anti-rabbit antibody, IgG-FITC (1:100, 
diluted with PBS (196 BSA), inspired with 488 nm waves), in 
the dark for 3 hours at 37.5 °C; rinsed with PBS (0.1% PVA) 3 
times and then stained in propidium iodide (Pl, 10 wug/mL, 
inspired with 568 nm waves) in the dark for 10 min; and finally 
rinsed with PBS (0.1% PVA) and fixed on slides for observing 
and photographing under a laser scanning cofocal microscope. 
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Statistical analysis 

Fluorescence densities were analyzed using LAS AF Lite and 
Image Pro-Plus. Every sample group was repeated three 
times and at least 20 zygotes or embryos at different stages 
were examined. Corresponding fluorescence data were 
obtained by running Image Pro-Plus. Data from alcohol 
treatment groups and control groups were compared with t- 
tests (Kim et al, 2002). GraphPad Prism 5 was used to 
produce histograms and fluorograms. 


RESULTS 


H3K9 acetylation of MII stage ovum 

Indirect immunofluorescence detection of H3K9 acetylation in 
MII stage ovum showed that no immunofluorescence could be 
found in the alcohol treatment groups or controls. However, 
immunofluorescences were found by DNA PI staining (data not 
shown). 


Immunofluorescence of H3K9 acetylation in pre- 
implantation embryos at different developmental stages 
in the in vivo groups 

Indirect immunofluorescence detection of H3K9 acetylation in 
pre-implantation embryos at different developmental stages in 
the in vivo groups is shown in Figure 1. 
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Figure 1 Different stages of mouse pre-implantation embryos 
and H3K9 acetylation immunofluorescence in vivo 

1-cell, 2-cell, 4-cell, 8-cell, morula and blastocyst in vivo treatment and 
control groups collected separately. Embryos were immunostained with 
specific antibodies against resistance to H3K9 acetylation antibody (H3K9 
in histone); in red (middle), DNA; Merged (right) . Scale bar=20 um. 


Semi-quantitative analysis of H3K9 acetylation in pre- 
implantation embryos at different developmental stages 
in the in vivo groups 

The analysis of fluorescence data showed that at the morula 
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stage, H3K9 acetylation levels in the in vivo alcohol treatment 
group were higher than in the control group, whereas at the 
blastocyst stage H3K9 acetylation levels in the control group 
were higher than in the alcohol treatment group. No 
differences were found for other stages. An abrupt reduction 
in H3K9 acetylation only occurred in the alcohol treatment 
group at the blastocyst stage. For other stages in both groups, 
H3K9 acetylation levels all initially decreased, then increased. 
High expression was found in the prokaryotic phase. 
Expression levels were rather low in the 2-cell stage and later 
increased. 


Control 
4 Treatment group 


Relative H3K9 acetylation 
fluorescence (100%) 





Different development stages in vivo embryos 


Figure 2 Semi-quantification of H3K9 acetylation levels of 
mouse pre-implantation embryos in vivo 


H3K9 acetylation intensities of mouse pre-implantation 
embryos in vivo from control group and treatment group 
quantified with Image Pro-Plus 6.0. Different columns and 
bars represent mean+SE, *: P— 0.05, control group in vivo 
compared with treatment group in vivo. 


Immunofluorescence of H3K9 acetylation in pre- 
implantation embryos at different developmental stages 
in the in vitro groups 

Indirect immunofluorescence detection of H3K9 acetylation in 
pre-implantation embryos at different developmental stages in 
the in vitro groups is shown in Figure 3. 


Semi-quantitative analysis of H3K9 acetylation in pre- 
implantation embryos at different developmental stages 
in the in vitro groups 

The analysis of fluorescence data showed that H3K9 
acetylation levels of the in vitro alcoholic treatment and control 
groups were comparable (except at the 2-cell stage) and that 
H3KO9 acetylation levels increased (Figure 4). 

Patterns of H3K9 acetylation in the in vivo and in vitro 
alcohol treatment groups were different (Figure 2, 4). At the 2- 
cell stage in vivo a reduction in H3K9 acetylation was found, 
and then increased. /n vitro, the H3K9 acetylation level 
increased, after which it decreased and increased again. For 
control groups, H3K9 acetylation levels in vivo at the 2-cell 
and blastocyst stages were lower than in vitro, while at the 4- 
cell-stage in vivo were higher than in vitro; no differences 
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Figure 3 Different stages of mouse pre-implantation embryos 
and H3K9 acetylation immunofluorescence figure in vitro 

1-cell, 2-cell, 4-cell, 8-cell, morula and blastocyst in the in vitro treatment 
and control groups were collected separately. Embryos were 
immunostained with specific antibodies against resistance to H3K9 
acetylation antibody (H3K9 in histone); in red (middle), DNA ; Merged 
(right). Scale bar-20um . 
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Different development stages in vitro embryos 


Figure 4 Semi-quantification of H3K9 acetylation levels of 
mouse pre-implantation embryos in vitro 

Different columns and bars represent mean+SE, control group compared 
with treatment group in vitro. 


were found for other stages. For alcohol treatment groups, 
H3K9 acetylation levels in vivo at the 1-cell and morula stages 
were higher than in vitro, while at the 2-cell and blastocyst 
stages in vivo were lower than in vitro. No differences were 
found for other stages. 

We found that the blastocyst rates of in vitro alcohol 
treatment group and their controls were comparable. 
However, for the morula and blastocyst rates of in vivo 
groups, the alcohol treatment group was lower than the 
control group. 


DISCUSSION 


Histone acetylation influences chromatin remodeling and 
gene expression. As one of the signs of gene transcription 
initiation, H3K9 acetylation functions as an important form of 
histone acetylation and changes in its level can directly affect 
gene expression, induce protein misexpression and lead to 
phenotypic defects or maldevelopment. Here, no H3K9 
acetylation indirect immunofluorescences in MII stage ovum 
were observed in either alcohol treatment groups or their 
controls and we assume this is because MII stage ovum were 
already in maturation and H3K9 deacetylated. 

In vivo results showed that H3K9 acetylation levels in the 2- 
cell stage were low, and then increased. The mechanisms 
responsible is likely that during the 2-cell stage zygotic genes 
are activated, accumulated activation effects suppress gene 
transcription, and then under the regulation of zygotic genes 
many development-related genes start to transcript (Wang et 
al, 2013). Pan et al (2013) claim that increased H3K9 
acetylation in mice embryos lays the ground for transcription 
initiation of development-related genes, which is in 
accordance with our data. From zygotes to the 8-cell stage in 
vivo, although H3K9 acetylation in the alcohol group was only 
slightly higher than the control group, a difference induced by 
alcohol occurred at the morula stage when many 
development-related genes would be activated. Under these 
circumstances, alcohol treated morula turned mature in 
advance and then cells began to deacetylate, and thereafter 
H3K9 acetylation levels were lower compared to the control 
group. Alcohol can harm female mice and affect H3K9 
acetylation of pre-implantation embryos in multiple ways. 
Studies on hepatic cells show that alcohol related high histone 
acetylation may be correlated with the MAPK signal pathway 
(Yeh et al, 2008). 

In vitro abnormal 4-cell stage H3K9 acetylation could be 
the result of activation delay of zygotic genes (Ma et al, 2001; 
Qiu et al, 2003). Zygotic gene activation (ZGA) is the initiation 
of embryo development from maternal regulation to embry- 
onic regulation when the embryo will start to synthesize its 
own mRNAs and proteins and no long depends on 
maternally-sourced mRNAs and proteins (Kidder, 1993). The 
2-cell stage is a vital period in mice embryonic development 
and is also the initiation phase of ZGA. /n vitro culture may 
delay ZGA and prolong the procedure from maternal 
regulation to embryonic regulation. 

Fluctuations in H3K9 acetylation levels observed in the vivo 
and in vitro control and alcohol groups from zygotes to 
blastocytes indicate that environmental differences can induce 
H3K9 acetylation imbalance in pre-implantation embryos. This 
imbalance may be correlated with the activation delay of 
zygotic genes and the postponing of 2-cell stage development; 
however, the mechanism for in vitro 2-cell stage development 
postponement remain unclear. Some studies have shown that 
at the 2-cell stage in vitro, the active oxygen content in 
embryos was much higher than that of other stages and the 
lesion caused by excessive active oxygen may be one of the 
reasons for 2-cell stage development postponement in mice 
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(Tang, 2006). Zhao et al (2005) reported that in vitro culture 
could suppress or decrease GCN5 and HDAC1 protein 
expressions in mice pre-implantation embryos and influence 
in vitro histone acetylation. Both in vivo and in vitro data 
suggest that long-term excessive alcohol consumption, even 
without alcohol intake during pregnancy, could still influence 
H3K9 acetylation in pre-implantation embryos because of 
remnant alcohol inside the body or metabolism products. 
Changes in H3K9 acetylation may also interfere with the 
correct establishment of epigenetic patterns in pre- 
implantation embryos. Alcohol may influence  H3K9 
acetylation in pre-implantation embryos and then induce gene 
misexpression by adjusting the activities of HATs; however, 
the underlying mechanisms at play here require further 
exploration. 
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ABSTRACT 


Since 1980, the white-browed crake (Porzana 
cinerea) has been experiencing an expansion from 
south of the Isthmus of Kra, northward to China. 
Recently, this species was observed in several 
locations throughout Southwest China, including 
Ningming and Baise, Guangxi (2012, 2013), and 
Xichang, Sichuan (2013). These sightings are the 
first distribution record of this species in mainland 
China, suggesting that the white-browed crake is 
following a natural species dispersal northward into 
mainland China from Southeast Asia. 


Keywords: Porzana cinerea, New record; Mainland 
China 


INTRODUCTION 


The white-browed crake (Porzana cinerea) is a small, slim- 
body crake with striking diagnostic face pattern that typically 
resides in vegetated coastal and freshwater wetlands, ranging 
from Malaysia and Greater Sundas through Philippines, 
Sulawesi, Moluccas and Lesser Sundas to New Guinea and 
north Australia, and east through Micronesia and Melanesia to 
Polynesia (Fiji and Samoa) (Taylor, 1996). The species was not 
included in the early avian checklist of China (Cheng, 1987), 
though the first recorded observation was in Southern Taiwan 
between 1984 and 1985 (Lin, 1997). Later, one vagrant was 
observed in 1991 at Mai Po, Hong Kong (Kennerley, 1992). To 
date, this species has never been recorded in mainland China 
(Wang et al, 2006; Zheng, 2011). During a recent survey of 
birds in Guangxi and Sichuan, several recent sightings in the 
area suggest that a new distribution record of this species in 
mainland China is warranted. 


OBSERVATIONS 


On 28 August 2012, one adult white-browed crake was 
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observed at Honggu reservoir, Ningming Town, Southwestern 
of Guangxi, China. This sighting was treated as vagrant 
record, and no photographic evidence was obtained or 
published. Several months later on 2 February 2013, a sole 
crake was found and photographed in a suburban pond 
locating in Baise City, in Northwestern Guangxi Zhuang 
Autonomous Region. According to this species’ typical face 
pattern and behaviour, it was identified as an adult white- 
browed crake. Later, one juvenile was found at farmland on 
15 November (Figure 1) and another at Longjing reservoir on 
24 December 2013 in Baise City, Guangxi. In addition, one 
juvenile white-browed crake was found by Min Xue, a 
birdwatcher, in Xichang Wetland Park, Sichuan Province, on 
25 October 2013 (Table 1, Figure 2). 





Figure 1 A white-browed crake (right) in Baise, Guangxi, South- 
western China (Photo by Dong-Sheng GUO) 
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Table 1 Records of white-browed crake in Mainland China 





Date Age class Number (n) Site recorded Longitude Latitude Altitude (m) Habitat 

28 Aug. 2012 Adult 1 ae E107° 09' N22° 03' 164 Reservoir 

02 Feb. 2013 Adult 1 Baise, Guangxi E106° 36' N23° 52' 142 Pond 

25 Oct. 2013 Juvenile 1 Xichang, Sichuan E102° 18' N27° 48' 1511 Wetland park 
15 Nov. 2013 Juvenile 1 Baise, Guangxi E106° 36' N23° 52' 147 Farmland 

24 Dec. 2013 Juvenile 1 Baise, Guangxi E106° 36' N23° 52' 142 Reservoir 
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Figure 2 The distribution range and records of white-browed crake 
m: New records for mainland China from this study; A: Previous records in Hong Kong and Taiwan; e: Records from non-Sundaic mainland Asia since 1980 
(Duckworth & Evans, 2007); grey shadow: The distribution range (Handbook of the Birds of the World, 1996); ?: Extinct site. 


DISCUSSION 


The white-browed crake was previously thought to distribute 
mainly in south of the Isthmus of Kra, which includes Malaysia, 
Singapore, Philippines, Indonesia, New Guinea and North 
Australia (Taylor, 1996). Accordingly, in Southeast Asia, early 
records were from Malay Peninsula (Robinson & Chasen, 
1936). In the 1990s, new distribution sites were found across 
Southeast Asia, including Thailand, Cambodia, Vietnam and 
Laos (Buckton & Safford, 2004; Mundkur et al, 1995; Robson, 
2000, 2004, 2011). Duckworth & Hedges (2007) analyzed this 
species distribution records since 1980 and proposed that 
white-browed crake has been expanding northward from a 
purely Sundaic distribution and moving almost to China. This 
hypothesis was well supported by the studied records: one 
distribution site in Laos was ~800 km north of any historical 
sighting, while the site in northern Vietnam at Van Long 
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Nature Reserve—quite close to the southwestern Chinese 
border—was considered as an “amazing range extension” 
(Duckworth & Hedges, 2007). Taking our distribution data of 
the crake in mainland China alongside several recent records 
for other birds in this region (Jiang & Ning, 2010; Jiang et al, 
2013), it seems likely that the  white-browed crake’s 
distribution range is quickly expanding in Asia. Indeed, the 
distance between Baise and Van Long Nature Reserve is 
~400 km, and Xichang to Van Long is ~900 km. 
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DEAR EDITOR: 


The  species-area relationship (SAR) describes the 
phenomenon whereby the number of plant and animal 
species found in an area of wild habitat is strongly correlated 
with the size of that area. As one of the few ecological laws, 
the SAR plays a vital role in the design and assessment of 
biodiversity protection regions (Lomolino et al, 2010; Ladle & 
Whittaker, 2011). 

Increasing area and habitat promote species richness 
(Triantis et al, 2003) and both these mechanisms have 
theoretical support. Basically, the notion of area per se claims 
that the effects of area are embodied in the chance of species 
extinction (Preston, 1960, 1962; MacArthur & Wilson, 1963, 
1967). In contrast, the habitat diversity hypothesis states that 
the highly diversified habitats often found in large patches 
positively impact species addition (Williams, 1964). Although it 
remains controversial which mechanism is more important in 
influencing species richness, some ecologists believe that the 
area per se and habitat diversity are mutually complementary 
and not mutually exclusive (Triantis et al, 2003). 

Land-bridge islands are characterized by well-delineated 
boundaries, an inhospitable surrounding matrix and relatively 
homogeneous habitats (Hu et al, 2011). Land-bridge islands 
are considered excellent models for studying the relationship 
between area per se and habitat diversity. Located in Chun'an 
County in Zhejiang, eastern China (N29°22'-29°50', E118?34"- 
119?15^, Thousand Island Lake is a large man-made lake 
following dam construction on the Xin'an River in 1959. With a 
water surface area of approximately 580 km, the lake 
contains 1 078 land-bridge islands larger than 0.25 ha (108 m 
in water level elevation). 

To understand correlations between area per se and habitat 
diversity we used bird data for 41 land-bridge islands in Thousand 
Island Lake from 2006 to 2009. This data was used to compare 
the goodness-of-fit of the classical Arrhenius SAR model, and the 
choros (K) model integrating both area and habitats (Triantis et al, 
2003). The canonical power function model (logS-logcx zxlogA, 
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in which S represents the number of species and A is the size of 
the area; c and z are constants) was used to investigate the effect 
of area on bird species richness. Choros (K) (K-HxA, in which H 
represents habitat diversity and A is the area of the study island) 
was adopted to assess correlations between area and habitat 
diversity (Zhang et al, 2008; Wang et al, 2010). The classical 
SAR model can be expressed as: 


Log (S)=log (c)+z log (A) (1) 
and the K model can be expressed as: 
Log (S)=log (c)+z log (K) (2) 


Akaike Information Criterion (AIC) was used to determine 
the optimum model and the model with the lowest AlCc 
(modification of AIC for small n) was considered better 
(Burnham & Anderson, 2002). Statistical differences in z 
values among regression equations were compared by 
referencing the methods of Zar (1996) f(b1—b2)/Sp1,52 in 
which bf and b2 are both regression coefficients and Sp1,52 
are standard errors of regression coefficients). 

Our results show that the lowest AlCc was found in the SAR 
model (Table 1), however, no differences in z values were 
found between the two models (P>0.05). These findings 
indicate that the Arrhenius SAR model gives a better fit than 
the K model and there is no effect of habitat diversity on the 
SAR. This finding is in contrast with that of Triantis et al (2003) 





Table1 | Comparison of the SAR and choros (K) models 
Models 2 Adjusted r2 AlCc 
(logS-logA) — 0.12 0.64 —219.71 
(logS-logK) — 0.11 0.62 —217.94 
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in which the K model better explained the effect of area and 
habitat diversity on species richness. 

One possible explanation for this discrepancy is that due to 
the relatively homogeneous habitats on the islands of 
Thousand Island Lake, insignificant differences in habitat 
diversity may be all that exists between the largest and 
smallest islands. Therefore, compared with changes in area, 
habitat homogeneity is more prominent in Thousand Island 
Lake and bird richness is strongly correlated with area (Hu et 
al, 2011; Ding et al, 2013). The weak effect of habitat diversity 
found in this study may also be the result of different statistical 
methods. There are various ways to define habitat diversity 
(Looijen, 1995, 1998); for example, Ricklefs (1979) defines 
habitat as vegetation coverage in a given area but Whittaker 
et al (1973) defines habitat as the multidimensional space 
occupied by plant and animal species. Although the latter 
definition (with minor improvements) has been widely 
accepted by ecologists (Krebs, 1988, 1994; Looijen, 1995, 
1998), the definition itself is too broad and has been difficult to 
apply consistently across different studies (Newmark, 1986). 

Here, only island area was important when explaining the 
relationship between area, habitat and species richness. 
Consequently, eliminating habitat diversity alters the results 
only subtly and in protection practice, area has the highest 
authority in maintaining bird richness. Effective evaluation of 
species biodiversity can be obtained by choosing area as a 
parameter, and priority protection areas can be determined 
accordingly. This finding suggests that more attention should 
be given to large islands with high species richness in 
Thousand Island Lake (Wang et al, 2010). 
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